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SYNOPSIS 


The work reported in the preeent thesis deals with the 
problems of flow and heat transfer through a porous medium. 

In recent years, considerable interest has been envinced 
towards the study of natural convective flow phenomena within a 
fluid-saturated porous medium; such studies are of intrinsic 
importance and relevance to the various applications in the areas 
of petroleum engineering, chemical engineering, biomechanics, and 
geophysics. 

For the modelling of flow through Isotropic and homogeneous 

\ 

porous materials, the basic constitutive relation which has been 
widely used so far is the Darcy's law which relates the seepage 
velocity to the pressure gradient. In many practical situations, 
however, the porous medium is bounded by impermeable walls and has 
a high permeability, which render the Darcy's equation inapplicable 
for flow modelling- The inclusion of viscous stress terms in the 
Darcy's equation giving rise to a boundary layer type of flow near 
the Impermeable wall, results in the so-called Brinkman model. A 
further extension of Darcy's law which takes into account the 
convective acceleration ' effects along with the viscous shear 

in the fluid is termed as the generalised Darcy's law. 

The oblectlve of the nresent work is to study two-dimensional 

using generalised Darcy's law as the 
idary layer situations as well as 
aerlcal techniques have been used to 


^n- linear governing equations. 



In the first chapter, an overview of the structure and 
properties of a porous medium is presented along with the 
conservation equations of mass, momentum and energy. Various flow 
models available for a fluid-saturat ed porous medium have been 
described. A brief survey of the published literature related to 
the present work is also presented. 

In the second chapter, the existence of similarity solutions 
for free convection along a vertical plate immersed in a porous 
medium is analysed- Necessary and sufficient conditions are 
obtained in a porous medium described by generalised Darcy’s law, 
the Brinkman model or the Darcy’s law. The approach consists of 
defining an independent variable and two dependent similarity 

functions corresponding to the velocity and temperature solutions 
in terms of generalised transf ormat ion functions. Substituting 
these generalised transformation functions into the boundary 
layer equations ^ a set of ordinary differential equations will 
result only if the coefficients of all the terms are made 

constants. The minimum set of conditions to be satisfied for 
making these coefficients constant are referred to as similarity 
conditions. For the particular situation with wall temperature 
varying as a power function of distance i numerical solution has 
been obtained using the shooting technique. Velocity and 

temperature profiles have been presented for different 
permeability parameters and Prandtl numbers. It is found that in 
the boundary layer^ the flow is weak and the temperature profile is 
less steep, when the voidage of the medium is small. Prandtl 

number is seen to have strong influence upon the temperature and 



velocity profiles due to its effect upon heat penetration into the 
.^Itiid. 

The third chapter deals with the problem of conjugate free 
convection between a vertical fin and a saturated porous medium. 
The governing equations are based on the generalised Darcy's model 
and are solved iteratively within the framework of boundary layer 
approximat ions by a highly-impl icit finite difference scheme. A 
variable mesh aiz^ is used in the direction normal to the fin 
surface in order to obtain the flow and heat transfer 
characteristics inside the boundary layer accurately. It is 

observed that the fin cooling is more effective at higher Grashof 
number or permeability parameter due to stronger convection 
effects. The local Nusselt number increases with the Grashof 
number, permeability parameter and Prandtl number and it decreases 
slightly with the conduct ion-convect ion parameter. 

Unsteady natural convection in a rectangular enclosure is 
presented in fourth chapter. The enclosure consists of two 
isothermal vertical walls and two adiabatic horizontal walls. A 
general parametric study in terms of the Rayleigh, Darcy and 
Prandtl numbers has been presented. . The flow and temperature 
field solutions are obtained using the finite difference 
technique. The upwind procedure is applied for convection terms 
along with the central differences for the diffusion terms for 
spatially discretizing the governing equations. The timewise 
discret izat ion is performed using a fully implicit scheme. The 
resulting equations are solved using 1 ine-^by-line technique 
coupled with the tridiagonal matrix algorithm. Streamlines and 
isotherms have been plotted for analysing the flow and heat 



transfer inside the enclosure at different times. The transient 


variation of average Nusselt number at the hot wall is also 
obtained for highlighting the rate of heat transfer between the 
vail and the fluid trapped inside the porous matrix. The flow is 
found to be unicellular and governed by the modified Rayleigh 
number; the product of the Rayleigh and Darcy numbers. 

The last chapter deals with the problem of free convection 
over a moving vertical plate in a porous medium. The problem has 
an exact solution for the values of Prandtl number greater than or 
equal to unity. However , for Prandtl number less than unity, there 
is no exact result. Therefore, for Prandtl number less than 
unity the problem is solved here analytically by the momentum 
integral method and numerically by finite difference method. The 
results obtained by the two methods agree well. The fluid 
velocity at a point decreases with Grashof number and Prandtl 


number . 
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CHAPTER I 


INTRODUCTION 

1.1 GEKIERXL BACKGROUND 

The sttidy of processes involving fluid flow and lieal transfer 
through saturated porous nedia has evoked inaense research 

interest in recent years. A broad range of applications including 
geothersal aystens, crude oil production, ground water pollution, 
fibre and granular insulation, filtration and storage of nuclear 
waste fflaterlal etc. have been modelled as porous media systems. 
In such processes the intersticial fluid trapped in the pores is 
subjected to vapourization, condensation and heat or mass 
transport due to the prevailing temperature or concentration 
gradients. 

Host o£ the analytical and numerical works in this area 
accept Darcy’s law [1] as the macroscopic equation of motion, 
which states that (seepage) velocity is proportional to the 
applied pressure gradient and that flow resistance is directly 
proportional to fluid viscosity and inversely proportional to the 
permeability of the me'dium. However, the Darcy's law becomes 
inapplicable when the porous medium has a high permeability and is 
bounded by impermeable walls. In these cases, due consideration 
needs to be given to factors such as inertial effects in the 
fluid. Later generalizations of the Darcy’s law have focussed 
attention on some of these factors. 

The inclusion of viscous shear stress terms in the flow 
equations and the consideration of a velocity field which 



satisfies no slip at iaperaeable valla, have been proposed by 
Brinkman [2]. This is the so-called Brinkman's extension, vhieh 
accounts for the flow resistance due to the fluid viscosity 
itself. Yamamoto and Ivamura [3] proposed a flow model in which 
the effects of inertia are taken into account via convection terms 
in the momentum equations. Obviously this is an important 
correction which has to be considered at higher flow rates through 
the porous medium. Yet another approach to handle the inertial 
effects was proposed by Forchheimer [4], which approximately 
accounts for the losses in the pressure due to flow separation in 
the pores, by the inclusion of an additional term in the flow 
equations. Uith all these extensions, the flow models for porous 
media situations have more or less the same general structure as 
that of Navier-Stokes equations for free fluid flows. 

Host of the studies available at present on porous medium 
applications do not consider a generalized Darcy's law. In the 
present work, this generalized Darcy's model is taken as the 
equation of motion to study two-dimensional natural convection 
problems with and without the use of boundary layer 
approximations. The governing partial differential equations for 
these problems are non-linear, which are amenable to analytical 
study only in some limited cases. In the present work, full 
numerical solutions have been obtained using finite differences or 
shooting techniques, along with analytical solutions where 
available. Detailed velocity and temperature distributions and 
Nusselt numbers signifying the overall heat transfer rates, have 
been presented. Comparisons with previous works have also been 
provided where possible. 



l.H DCSCRIPTIOM OF THE STRUCTURE OF A POROUS MEDIUM 


Exaatplea of porous aaterlals are nuaerous. A few coaaonly 
encountered porous aaterials are soil, porous or fissured rocks, 
ceraaics, filter paper and a loaf of bread. Somewhat leas obvious 
examples, but still a part of this group, are large geological 
formations where open passages aay be of substantial si; 2 e and far 
apart. All of these examples have some properties in common that 
permit them to be grouped together as porous media. 

A porous medium may loosely be described as a 'solid with 
holes’. However, a solid block with isolated holes or pores would 
not normally be classified as a porous medium for modelling flow 
through the block. However, the same block may be acceptable as a 
porous solid in connection with thermal insulation. In flow 
applications, the definitions may be somewhat improved by 
stipulating that the pores be interconnected by several continuous 
paths from one side of the medium to the other with an even 
distribution (in either a regular or random manner) of holes and 
paths over the entire porous medium domain. 

Summarising these preliminary observations, the porous medium 
is defined [5] as: 

(i) A portion of space occupied by heterogeneous or multiphase 
matter. It should include a solid phase and atleast another phase 
which is not a solid but may be a gas or liquid. The solid phase 
is called the solid matrix. That apace within the porous medium 
domain which is not a part of the solid matrix, is referred to as 
the void space (or pore space). 

(ii) The void space should be evenly distributed throughout the 
porous medium and should be present inside each representative 
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elementary volume. An essential characteristic of the porous 
medium is that the specific surface of the solid matrix (defined 
as the total intersticial surface area of pores per unit bulk 
volume of the porous medium) is relatively high. In many 
respects, this characteristic determines the interaction between 
the solid matrix and the fluid in porous media. Another basic 
feature of a porous medium is that the pores comprising the void 
space are relatively narrow. 

(iii) At least some of the pores should be interconnected. The 
interconnected pore space is sometimes referred to as the 
effective pore space. As far as the flow through porous medium is 
concerned, the unconnected pores may be considered as part of the 
solid matrix. Certain portions of the interconnected pore space 
may, in fact, also be ineffective with regard to the flow through 
the medium. For example, pores may be dead-ended so that no flow 
occurs through them. Another way to define this pore 
interconnection characteristic of a porous medium is by requiring 
that any two points within the effective pore space be connected 
by many curves with an arbitrary maximal distance between them. 
For a finite porous medium domain, this maximal distance is a 
macroscopic property which is dictated by the dimensions of the 
domain. 

1.3 KIFINITION OF P<»OUS MEDIUM PRC»»ERTIES 
1 . 3. a Poromlly 

Porosity is defined as the fraction of the bulk volume of 
material which is occupied by voids. Thus, the porosity ^ is 


given by 



( 1 . 1 ) 


V 

_ P _ Volune o£ pores 

^ Bulk volume 

Since the portion of bulk volune not occupied by porea ia 
occupied by aolid graina or the natrix of the naterial, it followa 
that , 

V 

_ _ S _ Volume of aolida 

^ ~ * ~ Bulk volume ' (1-2) 

From the above definitiona, it ia evident that ^ is a macroacopic 
property. Two claaaea of poroaity can be defined, namely, the 
abaolute or total poroaity and the effective porosity. Absolute 
porosity ia the fractional void apace with reapect to the bulk 
volume, regardless of the pore connections. Effective poroaity is 
the fraction of bulk volume constituted by interconnecting pores. 
Hany naturally occurring rocks such as lava and other igneous 
rocks, have a high total porosity but essentially zero effective 
poroaity. Effective poroaity la an Indication of permeability of 
the porous medium but not a measure of it [6]. 

1 . 3. b Pmrmmabl lity 

Permeability is that property of porous material which 
characterises the ease with which the fluid may be made to flow 
through the material by an applied pressure gradient. 
Permeability can also be thought of as fluid conductivity of the 
porous material. 

The value of permeability is determined by the structure of 
the porous material. It has the dimensions of the square of 
length and is roughly a measure of the mean-square of the pore 
diameter in the material. Many porous materials have a 
directional quality in their structure. As a consequence, the 
permeabilities measured with flow perpendicular to each face of a 
cube may not be equal. Such materials are termed as anisotropic. 



1.4, GOVERNING EQUATIONS 


The analysis of flov and heat transfer In a porous medium is 
based on transport equations resulting from the differential 
balance laws. In many situations, detailed information of the 
velocity and temperature field is required to predict global 
effects such as flow resistance or heat flux from a given object. 
The flow and temperature fields, in turn, are extracted from the 
solution of the associated transport equations subject to 
pertinent boundary conditions. Uhen flow occurs through a complex 
structure such as a porous medium, the transport equations are 
still valid inside the pores, but geometric complexity prevents 
general solutions of the detailed velocity and temperature fields. 
Since it is impossible to know what is happening in each of the 
many pores, at best one can attempt to predict the statistical 
average of the physical quantities over certain representative 
volumes of the system. This aspect has been investigated in 
detail by Slattery [7,8] and Uhi taker [9]. They have considered a 
representative volume that is large compared to the pore size of 
the materal but small compared with the characteristic length of 
the object. 

Here, we shall define the characteristics of an idealized 
porous medium which is often employed in the modelling of porous 

media flows. It is assumed that the medium is isotropic and 

consists of fine particles which are uniformly scattered and fixed 
Jn space. Shrinking and swelling of the porous medium during the 
flow is neglected and the density of the particles is assumed to 
be constant. Any suitable microscopic volume of the medium 

contains both solid and fluid phases. In other words, the matter 



in any J nflnitesioal volune of. the ayetea under consideration will 
always be a solid-fluid mixture. For the fluid phase, a slnftle 
liquid or eas Is considered. Thus the state of an infinitesimal 
volume of the phases will be defined by the properties of the 
individual phases for the sake of simplicity in modelling. These 


properties are called intrinsic 

or 

actual properties. 

It 

is 

soiRetliaes assuiiied that 

each of 

the 

pores occupies the 

entire 

control voltiflie at any 

given 

point 

in the system. 

It 

is. 

therefore, necessary to 

modify 

some 

of the properties 

of 

the 


phases by appropriate volume averaging. These modified properties 
are called the bulk properties. From the above considerations, it 
is evident that for applying thermodynamic and hydrodynamic 
equations to an infinitesimal volume of the system, it will be 
necessary to use bulk properties of the phases [10 J. 

1.-4. a Consarvablon of Mass 

As discussed above, an infinitesimal volume of the system at 
any point is assumed to contain both the phases and each phase is 
assumed to occupy the entire volume. The velocity of the fluid 
phase is therefore defined as the fluid mass flow rate per unit 
area of the system (usually called the mean filter velocity) at a 
given point divided by the bulk density of the fluid at that 
point. Since a continuum model is considered, in general, for a 
three-dimensional averaged flow the mass conservation statement 
reads as : 

|| + = 0, (1-3) 

where, y la the volume averaged velocity vector. For an 
incompressible fluid, one may write. 


V.v 


0 


(1.4) 



1 . 4. b ConsArvA-tlon of" MoMonlink 


Darcy * a M odal 

Moat o£ the available analytical work on fluid flow and heat 
tranafer through poroua media are baaed on Darcy'a law. Thia law 
waa eatabliahed in 1856 by Darcy empirically, and aince then haa 
been verified experimentally by numeroua inveat igatora for 
Newtonian fluid motion in poroua media at low flow ratea. It 
atatea that for an iaothermal fluid moving with a alow ateady 
velocity u under a pressure gradient ^ through a poroua 
homogeneous and iaotropic bed. 


u 


K ^ 

/j dx ’ 


Cl-5) 


where, is the fluid viscosity and K is an empirical constant 
called the permeability of the medium. Originally the Darcy's law 
waa established to obtain the average or apparent flow rates of 
fluids through porous media in one-dimensional flow. Later 
researchers, assuming the fluid velocity to be a continuous 
function of space co-ordinates have given Darcy's law for flow 
fields in two or three dimensions. The extension of Darcy’s law to 
multi-dimensional flow including the body force term is expressed 
as. 


[-'TP * t} , 


( 1 . 6 ) 


where, t is the body force vector given by the weight of the fluid 
per unit volume. It is to be noted, however, that Darcy'a law la 
established empirically and it can not be derived analytically, by 
performing a momentum balance on a small element of the porous 
medium. It is merely a constitutive relation which does not yield 
much information about the permeability of the poroua medium. The 


drawbacks of Darcy's law are that it does not satisfy the no slip 
condition and is not valid for high flow rates. 

Non-Darcy Models 

There is enough evidence to indicate that Darcy's law is 


valid as long as Reynolds number (R - 


uK 


1/2 


) is small compared to 


unity. The frictional resistance offered to flow by the solid 
matrix, is the predominant factor in this flow regime. As the 
velocity of the flow increases, a region of gradual transition la 
observed from laminar flow with dominant viscous forces, to still 
a laminar flow with inertial forces governing the flow. Often the 
value of Reynolds number R = 100 is mentioned as the upper limit 
of this transition region in which Darcy's law completely breaks 
down. Some researchers explain this deviation from the linear law 
by the separation of flow from the walls of the solid matrix 
caused at large Reynolds number by inertial effects. 


( i ) Forchheimer ' s Model 

There Is no universally accepted non-linear law which is 
valid for Reynolds number greater than unity. Many such 
relationships appear in the literature. The more general law 
follows from Forchheimer ’ s modification [11] of the Darcy's law 
which may be written as 

= ^- ■‘^^'^1-1- (1-7) 

where, b la another empirical constant. Ergun [12) and Hard [13] 
suggested a similar equation of motion. 

( i 1 ) Br 1 nkman Model 

In Darcy’s law the term representing the viscous stresses In 
the fluid is not considered. This may lead to significant errors 


near an iaparaeabla solid boundary vhsre a boundary layar with 
large velocity gradients may exist due to the no slip condition at 


the wall. Thus, for 

accurate modelling 

one 

has 

to 

take 

into 

account the 

vlacoua 

ahear 

stresses, in 

addition 

to 

the 

Darcy 

reaiatance , 

however 

aiiiall 

they may be in 

the 

bulk 

o£ 

the 

fluid 


flow. Such an approach was postulated for the first time by 
Brinkman [2], for which a rigorous theoretical Justification was 
given later by Tam [14] and Lundgren [15], Their contention was 
that at the rigid boundary the no slip condition must be applied 
and this is inconsistent with the use of the usual Darcy equation. 
Based on these considerations, the momentum equation proposed by 

Brinkman is of the form: 

^ 2 

VP = - — v + #jVv, (1.8) 

where, fj is the viscosity of .the fluid. 

( i J i ) Generalised Darcy's Law 

In situations such as the flow through porous filters used 
for separating suspended particles from a gas stream, the inertial 
terms are quite Important. Several studies [16-20] have been made 
to analyse such flows on the basis of a generalised Darcy’s law 
which bears a lot of similarity to the Navi er-Stokes equations for 
the free fluid. 

Nost porous filters are made of fibrous materials and their 
porosity Is close to unity [21,22]. In order to accurately 
predict the flow In such a highly porous medium under 

non-lsothermal conditions, the buoyancy effects needs to be taken 
into consideration. Also, while accounting for the viscous 

stresses in the fluid due to velocity gradients in the flow, the 

form of the viscous term may be taken as in the case of a pure 



f3uld. Thus, the fundamental equations for an incompressible flow 
through a medium of high porosity are given by: 

[ yf + (y.V) v] = - VP+pi7v + pg-^v. (1.9) 

Uhen Reynolds number R << I , the convection terms can be safely 
neglected and equation (1-9) reduces to the one proposed by 
Brinkman. Also if the porosity of the medium is low, the Darcy 
term will dominate and it will chiefly account for the pressure 
drop. Thus, equation (1.9) may be taken as the generalized 
momentum equation which boils down to the equations corresponding 
to the other models in the appropriate range of parameters. 

1.4. c Cons*rvat.lon of' Enorgy 

The energy equation for a fluid-saturated porous medium was 
obtained by Uooding in 1957 .[23]. For an isotropic homogeneous 
porous medium in three-dimensional flow, this equation can be 
written as: 

H + (v.V) T - « (1.10) 
where, T is the temperature and « la the thermal diffusivlty of 
the fluid-saturated porous medium. In this equation, the viscous 
dissipation effect and the internal heat generation have been 
neglected. It is to be noted, however, that the above form of 
energy equation for a porous medium is derived by averaging the 
energy equations for fluid and solid parts over an elementary 
voi ume . 

l.S LITERATURE SURVEY AM) THE PRESENT 1PORK 

Early theoretical and experimental studies of convective flow 
of a fluid in a porous medium were reported In connection with the 
salt in subterranean sand layers [ 24-27 ] . Study of the motion of 
underground water with particular reference to geothermal 



phenomena was presented by Uoodlnp. [23,28-33], Eldet [34-37] and 
McNabb [38]. The results were presented on the structure of flow 
field and the correspondinft heat transfer rates fof convection of 
fluid in a porous medium owing to heat generatio*' from below. The 
criterion for the onset of convective flow predicted 

theoretically by Lapwood [39], while Uooding [33] McNabb [38] 

introduced the concept of boundary layers in porous media. 

Boundary layer analysis has been applied to problem of 

free convection in a porous medium about a vertical flat plate, 
when the wall temperature Is taken to be a pousf function of 
distance from the leading edge, by Cheng and They 

obtained similarity solutions for the problem. ^ subsequent 

paper Cheng [41] used the Karman-Pohlhausen integral method to 
obtain the local Nusselt number for a number of convective heat 
transfer problems in a porous medium. It was found that the 
approximate solutions for the Nusselt number obtained from the 
Integral method with different assumed profile* are in good 
agreement with those obtained by similarity solutions. The 
problem of natural convection adjacent to a vertical or horizontal 
heated surface with a power law variation in the wall temperature 
was studied by Cheng and Chang [42]. The results were obtained by 
the perturbation method where the first approxi®®^ gives the 
boundary layer solutions. A numerical solution for free 

convective flow past a vertical semi-infinite fist plate embedded 
in a highly saturated porous medium with a non-uniform temperature 
or a non-uniform heat flux distribution on surface was 

developed by Na [43]. Both local heat transfer rate and excess 
surface temperature as a function of distance along the plate were 
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tabulated for a few oaaea of preacrJbed wall tenperature and heat 
flux d latribut iona . A atudy of the flow and heat tranafer along a 
vertical flat plate which la placed in a aaturated poroua aediun 
and at an arbitrary diatance above the horizontal wall waa 
performed by Ingham and Pop [44], The firat and aecond order 
boundary layer equationa and the correapondlng outer Invicid flow 
equationa were atudied to find the effecta of large, but finite 
valuea of Rayleigh numbers. Results were obtained when the plate 
la maintained at either constant temperature or uniform heat flux. 
An approximate aeries solution was found for free convection 
boundary layer flow past a vertical plate in a porous medium, over 
a wide range of wall temperature and Injection velocity 
conditions, by Govlndara julu and HolarvJzthi [45]. Although the 
treatment waa similar to that of Herkin [46] the governing 
equationa and boundary conditions were different. 

The problem of conjugate heat transfer from a downward 
projecting fin was studied by Pop et al . [47, 4&]. In such 
problems, inclusion of the coupling between the surrounding fluid 
and the fin makes the analysis very different from that of a 
vertical plate. Similarity solutions of the appropriate boundary 
layer equations were presented . 

The buoyancy induced boundary layer flow adjacent to 
horizontal impermeable surfaces in a saturated porous medium was 
performed by Cheng and Chang [49] . Similarity solutions were 
obtained for the convective flow above a heated surface or below a 
cooled surface, where the wall temperature is a power function of 
distance from leading edge. Applications of the problem to 
convective flow in a liquid dominated geothermal reservoir were 



Steady free 


convection in a porous 


also discussed. Steady free convection in a porous medium 

adjacent to a horizontal impermeable surface with prescribed wall 
temperature distribution was investigated by Chang and Cheng [SO], 
using the method of matched asymptotic expansions. Similarity 

solutions for buoyancy induced flow in a saturated porous medium 
adjacent to a horizontal impermeable surface were given by 
Chandrasekhara {51]- The analysis incorporates the variation of 
permeability from the wall. An analytical study was presented by 
Ingham and Pop [52] for moderately high Darcy-Rayleigh numbers 
for steady free convection adjacent to a horizontal upward- facing, 
heated surface which is bounded by a vertical wall. The solution 
of the governing equations was obtained by the method of matched 
asymptotic expansions. 

The transient free convection about a vertical flat plate 
embedded in a porous medium adjacent to a semi-infinite flat pi 
with a step increase in the wall temperature and the surface heat 
flux was initiated by Cheng and Pop [53]. The governing equations 
are first order partial differential equations of hyperbolic type 
and these have been solved exactly by the method of integral 
relations. Similarity conditions and similarity equations for 
steady and unsteady natural convection adjacent to a semi-infinite 
inclined flat plate in a porous medium were examined by Johnson 
and Cheng [54]. Similarity solutions were obtained by Cheng [55] 
for mixed convection over a horizontal impermeable surface in 
porous medium. The problems of free and forced convection 

boundary layer flows in porous media, have been dealt by Cheng 
[56] on the basis of boundary layer approximations. Similarity 

solutions were obtained for the special case when the far 



velocity and the wall teaperature dJ atrJ but 1 one vary according to 
the aaae power function of distance. The free convection boundary 
layer on a aenl-lnf inite heated flat plate which Ja inclined at a 
amall angle to the horizontal waa considered using the aeries 
aethod by Inghan et al . 

The analysis dealing with the onset of convection in a fluid 
saturated porous box, heated from below was performed by Lowell 
[58]. In an experimental study of natural convection in 
liquid-saturated confined porous medium, Kaneko et al . [59] had 
shown that the mode and intensity of convective motion are 
affected by the angle of inclination of the medium and the 
properties of the fluid- The true critical Rayleigh number for the 
case of a convective flow of a fluid in a rectangular box of 
porous material heated from below waa found for various box 
geometries by Beck [60], The preferred cellular mode of the 
motion at Rayleigh number Just above the critical value was 
determined and the invalidity of Darcy's law in determining the 
stability of flow was also discussed. Natural convection in a 
confined porous medium was investigated experimentally and 
numerically by Holst and Aziz [61, 62]. 

Buoyancy-driven convection inside a porous cavity due to 
heating in the horizontal direction was analysed by Ualker and 
Homsy [63]. The solutions were obtained over a range of the 
dimension] ess parameters of Darcy-Rayl elgh number and the cavity 
aspect ratio. Finite amplitude cellular convection in cubic boxes 
containing fluid- saturated porous material were reported by 
Strauss [64] for Rayleigh number as large as 150. Steady natural 
convection in a slightly inclined rectangular box filled with a 



ID 


porous asdiuffi v&a studied by U&lch and Dulieu [65] using the 
Galerkin technique. Exact aolutiona were developed by Philip [66] 
for anal 1 cavities of various shapes due to a uniform temperature 
gradient normal to the gravitational field. Exact solutions were 
presented in detail for rectangular and elliptical cavities for 
arbitrary aspect ratios and also for equilateral triangular 
cavities. Simple exact solutions have been found to exist for 
cavities of other shapes also. The finite difference technique 
was used to investigate finite amplitude convection in a porous 
container with a fault - like geometry. The principal purpose was 
to determine the role of boundary conditions and the initial 
perturbation to the subsequent flow pattern. The results were 
compared with data from naturally occuring geothermal systems by 
Lowell [67]. The thermally driven flows inside a two-dimensional 
rectangular cavity filled with a fluid saturated porous medium was 
considered in the large Rayleigh number limit, when the applied 
temperature gradient is perpendicular to the gravity vector. 
Bounds on the extent of various flow regimes In the parameter 
space were deduced by Blythe et al.[6S]. A numerical and 

experimental study of two-dimensional steady free convection was 
presented for various thermal conditions on the walls of a 
rectangular cavity by Prasad et al-[ 69-71]. 

Buoyancy-driven convection in a differentially-heated 
vertical porous layer was studied theoretically by Gill [72] and 
Weber [73]. The model was of finite extent and the temperature 
difference between the vertical walls was assumed to be large. 
Satisfactory agreement with experiment was obtained for the 
interior temperature distribution and the Nusselt number. Gill 



dev«3op«d & modified Oaeen technique and Ueber, in his atudy, 
extended the method to Include some effecta of kinematic 
viacoalty. Bejan [74] altered GJll'a theory by calculating the 
net flow of energy (conduction and convection) in the vertical 
direction and aetting it equal to zero at the top and bottom enda. 
In a further atudy Bejan [75] modified Ueber'a theory for 
convection in a vertical poroua layer by combining Ueber'a 
boundary layer aolutlon with zero energy flux boundary condltiona 
at the top and bottom aurfacea of the alab. Reaulta were alao 
obtained by Bejan and Tien [76] to conalder the corner effecta 
particularly in aha! low cavitiea. Several vertical layera with 
different permeabll it lea were atudied by Poulikakoa and Bejan 
[77]. A fundamental atudy of natural convection in a fluid 
aaturated poroua layer differentially-heated in the horizontal 
direction, when the fluid exhlbita a denaity maximum waa preaented 
by Poulikakoa [78]. Reaulta were obtained in [79-82] for the 
varioua poroua layer problems with different condltiona. 
Convection in a tilted layer waa atudied by Ueber [B'^J while 
Blythe and Simpkina [84] atudied the convection for a temperature 
dependent viacoalty. 

Forchheimer extenaion of Darcy'a equation of motion waa first 
used by Plumb and Huenefeld [85] to study free convection from a 
vertical plate embedded in a porous medium. Be jan and Poul ikakoa 
[86] reconsidered the problem and preaented boundary layer 
solutions for the intermediate and non-Darcy flow regimes to 
supplement the boundary layer analysis of Cheng and Mlnkowycz [40] 
for Darcy flow. Their reaulta indicate that the slope of the 
curve of Nuaaelt number versus Rayleigh number decreases when the 



flow regime changea from Darcy to non-Darcy behaviour . 
Poulikakoa and Bejan [S7], Potillkakoe [BS] and Praaad and Ttintomo 
f893 also used the Forchheimer model to study non-Darcy natural 
convection in a vertical cavity for isothermal and constant heat 
flux conditions on the side walls. They have shown that, a new 
dimensionless parameter called the Forchheimer number governs the 
non-Darcy flow through porous media and the overall Nusaelt number 
decreases with an increase in Forchheimer number. 

It is noteworthy that with a few exceptions, most of the 
published theoretical work on convective heat transfer use Darcy's 
law as the momentum equation. Mathematically, since the order of 
Darcy's law is lower than the Navier-Stokes equations, the no slip 
boundary condition can not to be imposed on the impermeable 
boundary. For flow through a porous medium with a high 
permeability , the Brinkman model was used by Chan et al . [90] to 

study convective heat transfer in an enclosed, gas-filled, porous 
medium. It was also used by Katto and Masuoka [91] as well as by 
Ualker and Homsy [63], to study the onset of convection of a 
liquid in a porous medium bounded between parallel plates and 
heated from below. The Brinkman model was employed by Hau and 
Cheng [92] to study the natural convection about a semi -inf ini te 
vertical plate In a porous medium. Flow of a viscous fluid at 
small Reynolds number past a porous body was investigated by 
Yamamoto [93]. Sen [94] studied natural convection in a porous 
cavity . 

Recently several works have been published using the 
generalised Darcy 's law, in which the convective acceleration and 
the viscous stresses are also taken Into account [19, 9S]. This 



general i7.ed Darcy's law was used by Rapt J a et al . [96] to analyse 

the steady free convection and heat transfer in a porous aedium 
bounded by a vertical infinite plate. The solution was obtained 
by the series aethod when the surface temperature has a constant 
suction velocity. Chawla and Singh [97] studied the oscillatory 
flow past a porous bed. Rapt is [98] analysed the unsteady free 
convective flow through a porous medium bounded by an infinite 
vertical plate, when the temperature of the plate la oscillating 
with time about a constant non-zero mean. Analytical expressions 
for the velocity and temperature fields are derived by using the 
series method. In a subsequent paper, Raptis and Singh [99] 
obtained numerical solutions of unsteady free convective flow 
along a moving vertical wall embedded in a porous medium. In a 
recent paper Raptis and Perdlkis [100] have studied combined free 
and forced convection through a porous medium. 

The present work concentrates on the analysis of 
two-dimensional free convective flows with and without boundary 
layer approximations. The generalized Darcy’s law is considered 
as the equation of motion. Since the governing equations are 
non-linear and there exists a mutual coupling between them, the 
solutions are obtained with the help of numerical techniques. 

The first problem considered here, deals with the possible 
similarity solutions for free convective flow on a semi-infinite 
vertical plate embedded in a saturated porous medium. The 
necessary and sufficient conditions are obtained for similarity to 
exist, based on a generalised similarity method. Numerical 
solutions using the shooting method are obtained for a particular 
case and the corresponding velocity and temperature distributions 



are preB<^ntftd. 

The third chapter diecnaaem the problem of conjugate free 
convection about a rectangular fin embedded In a poroua medium. 
The governing equations are solved by using a highly-implicit 
finite difference scheme with a variable mesh-size in the 
direction normal to the fin. Uall temperature dJ stribut ions and 
the local Nusselt numbers are obfiained and plotted for various 
parameters of the problem. 

A parametric study has been presented in chapter four for 
unsteady free convection inside a rectangular enclosure with two 
isothermal vertical and two adiabatic horizontal walls. The 
velocity and temperature distributions are obtained using the 
finite difference technique. The upwind-di f f erencing procedure is 
followed for convective terms to get the spatial discr etizat ion of 
the governing equations. A fully implicit scheme is used to get 
the timewise discretizat ion and the resulting algebraic equations 
are solved by the line-by-line procedure coupled with the 
tridiagonal matrix solver algorithm. 

The last chapter deals with free convection about a moving 
vertical wall embedded in a porous medium. The problem la solved 
analytically by the momentum integral method and numerically by 
the semi-implicit finite difference method . The results obtained 
by the two methods are in good agreement. The veloc ity 
diatribution and akin-friction at the wall are obtained. 
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CHAPTER II 


POSSIBLE SIMILARITY SOLUTIONS FOR FREE CONVECTIVE FLOW ALONG 

% 

A VERTICAL PLATE IMMERSED IN A POROUS MEDIUM 

2.1 IMTRODUCTICW 

There has recently been considerable interest in the study of 
free convective flows in porous media because of their 
applications to many physical and engineering problems. 
Theoretical studies in this area leading to exact solutions, 
chiefly pertain to the linear stability analysis of the onset of 
free convection [1-3] and the prediction of velocity and 
temperature distributions for free convective flows near heated 
surfaces [4,5]. 

The problem of natural convection in a free fluid has been 
investigated by many researchers and notably among them are 
Stewartson [6], Gill et al , [7], Rotem and Clausen [8], Pera and 
Gebhart [9] and Blanc and Gebhart [10]. In the above 

investigations, similarity solutions have been obtained using 
boundary layer approximations. For flow through a porous medium 
also, Uooding [11] and IfcNabb [12] showed that the convection 
takes place in a thin layer around the heated surface and the 
boundary layer approximations are valid. 

A large volume of literature exists on thermal boundary layer 
studies for natural convection inside a saturated porous medium 
using Darcy’s law as the momentum equation. Cheng and Hinkowycr. 
[13] have considered the idealized problem of a vertical flat 
plate embedded in an isotropic porous medium with the wall 



t eisperattire preacribed as a power function of heii^ht. Uithin the 
framework of boundary layer approximat ions the similarity 

solutions were obtained. In another paper Cheng and Chang [14] 
have obtained the similarity solution near a horizontal plate. 
The similarity solution for combined free and £orc€td convection 
about an inclined surface in a saturated porous medium were 
analysed by Cheng [15]. Considering the effects of lateral 

injection or withdrawal of fluid through the wall, Cheng [16] has 
analysed possible similarity solutions for power law type 
variation of wall temperature. In a subsequent publication [17], 
he has also found the similarity solution for mixed convection 
over a horizontal impermeable surface in a saturated porous 
medi um. 

Na and Pop [16] have investigated steady natural convection 
heat transfer between an impermeable vertical surface and a porous 
medium, when a non-uniform heat flux condition is prescribed at 
the wall- They have shown that the governing equations possess a 
similarity solution for only a limited class of wall temperature 
or heat flux distributions. Transient boundary layer flow in a 
porous medium is dealt by Johnson and Cheng [19]. Similarity 
solutions have been obtained for specific variations of wall 
temperature with time and position. Plumb and Huenefeld [20] and 
Be Jan and Poulikakos [21] conaidered a non-Darcy model to study 
free convection from heated surfaces in a porous medium. They have 
shown that only for the case of a vertical isothermal surface 
similarity exists. All of the above theoretical solutions based on 
Darcy or non-Darcy models, do not take into account viscous force 



acting along the impermeable surface and hence the no slip 
condition at the boundary is not satisfied by them. 

The objective of the present work is to obtain the necessary 
and sufficient conditions under which similarity is possible for 
free convection adjacent to a vertical plate. In a porous medium 
described by the generalized Darcy law. The generalized similarity 
method which was earlier used by Tang [22], Cheesewrlght [23], 
Gebhart [24] and Johnson and Cheng [19], has been adopted here. 
The approach consists of defining an independent similarity 
variable r) and two dependent similarity functions of 17 
corresponding to the velocity and temperature solutions, in terms 
of generalized transformation functions. Uhen these variables are 
substituted into the boundary layer equations they can be reduced 
to a set of ordinary differential equations only if the 
coefficients of all the terms are made constant. The minimum set 
of conditions to be satisfied for making these coefficients into 
constants are referred to as similarity conditions. In general, 
the similarity conditions contain derivatives of the 

transformation functions and are thus differential equations 
themselves. Any set of solutions to these differential equations 
which will result in explicit forms for the transformation 
functions, can be used to transform the boundary layer equations 
of the problem into a similarity form. 

The first part of the present chapter deals with the 

theoretical discussion and the derivation of the similarity 
equations. In the later part, the numerical solution of the 
situation with the wall temperature varying as a power function of 
distance, has been obtained. 



2.2 I30VERNING E<»lATIC»ftS 

To aimplify the math emat leal analyais, the following 

asstimpt ions are made: 

(1) The fluid and the porous medium are homogeneous and 
isotropic. 

(2) The porojua medium is fully saturated. 

(3) The equations of motion are based on generalized Darcy's 
model - 

(4) The Boussinesq approximation is valid, i.e. the fluid density 
can be taken as constant except in the buoyancy force term* 


(5) Other fluid and porous medium properties are constant. 

Uith the above assumptions the boundary layer equations for 
laminar free convection on a vertical plate are given as: 
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( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 


where, u and v are components of velocity in the tani^ential and 
normal directions to the plate, g is the acceleration dtie to 
gravity, ft is the coefficient of volume expansion, v is the 
kinematic viscosity, T is the temperature of the fluid, T is the 
ambient fluid temperature, is the permeability of the medium 
and cn denotes the thermal diffusivity. The variables x and y 
denote space co-ordinates in the horizontal and vertical 
directions. The physical geometry of the problem is given in 


£ igure 1 . 



The situation, considered for the present analysis 
corresponds to greater than T^. By a siiaple change in 
co-ordinate system, one can get the results for the case when T 

w 


is smaller than T also, from the present result. The initial, 
conditions, depend on the problem under considerat ion, while tihe* 
boundary conditions are as follows: 

y=0: u=v=0, T=T (x,t), (2.4) 

w 

y-*-ao: u = 0,T=T (constant), (2.5) 

and the boundary conditions at the leading edge are: 

X = 0 : u = V = 0, T = T (2.6) 

Ob 

Introduce the following dimensionless quantities, 
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X y uL vL ^ x>t 

i< 


g/?I, (T-T ) 

Ob 


K = — - , Pr = v/<*, 

I. 

(2.7) 

where, L is the characteristic length, Pr is the Prandtl mimber- 
and G is the generalised Grashof number. 

In terms of the dimensionless quantities of equation (2.7), 
the equations (2.1)-(2.3) reduce to the following non-dimenai onal 


( 2 . 8 ) 

(2.9) 

(2.30) 

( 2 . 11 ) 
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The boundary conditions transform to 

y=0 : u=v=0, G=G (x,t), 

w 



y->.<»:ij = 0,G = 0, (2 ,“1 2) 

X = 0 : u = 0, G = 0. (2.13) 

The continuity equation (2,8)i8 automatically satisfied when a 
dimensionless stream function is introduced such that. 


Sy 


V = 


Wx 


(2.14) 


Ue shall choose a similarity variable y) in the following form, 

n = y (x,t) (2.35) 

and the corresponding dependent variables i(n) and B(i^) as, 

f(r)) - ry,t) 6 ( 7 ?) - 

2 w 

The present problem of finding the conditions for the existence of 


si mi 3 ar ity , therefore , 


reduces to that of determining the 


functions ^,(X,t), ^-(x,t) and G (x,t) such that the equations 
1 Z V 

(2.9) and (2.10) are transformed into ordinary differential 
equations for £(7}) and ^( 17 ) with appropriate boundary conditions 
for these variables. Combining equations (2.15) and (2.16), we 
8et , 


and 


u = ^ ^ ^ f 

^3 ^2 

^ ^2 ^2 

V = - ^ = - r ( £ ) f + — ( - 

^ ^x 


= G e. 
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) 7? f’ 1 


(2.37) 

(2-18) 

(2-19) 


Substituting these expressions in equations (2.9) and (2.10), the 
following generalised siaailarity equations are obtained: 


£”’ - a^T? £” - (ajj+a^+ag) £’ + a^ £ f 


(a +a ) £’ + a d = 0 

4 5 1 


( 2 . 20 ) 


JJ, 8” - ca^n - a^f) 8- - (a^ * a,f-) 8 = I) 
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( 2 . 21 ) 



where, priwe denotes derivative with respect to 17 . 
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( 2 . 22 ) 

(2.23) 

(2.24) 

(2.25) 

(2.26) 

(2.27) 

(2.28) 

(2.29) 


It may be noted from equation (2.29) that the permeability 
parameter, K may be a function of space and time both. The 
boundary conditions for the variables f and B are: 

f(0) = £’(0) = 0 , £’(<») = 0 (2.30) 

d(0) = 1 , d(a>) = 0 . (2.31) 

It is clear from equations (2.20) and (2.21) that the 
similarity solutions are possible if and only if a^ , a^ , a^, a^ , 

a^ , a,, a_ and a^ are constants. Consequently, equations ( 2 . 22 ) - 

bo/ o 

(2.29) may be regarded as the necessary and sufficient conditions 
for similarity solutions. Since there are more conditions than 
unknovns , the constants are not all arbitrary. For example, the 



following relations can be easily verified 


= ^^2 ^3 

= 3a, + a,. 
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(2-32) 

(2.33) 


A general method for the determination of the unknown functions 
and G ,is to solve for and 4^ first from the equations 

1 4* V X 4U 

(2-23) - (2-26), and then solve for the surface temperature 
variation from equation (2.22). 

2.3 SIMILARITY SOLUTICM4S 

Case - 1: Steady free convection with constant surface 

temperature . 

By 


For this case, a^ = a^ = a^ = a^ 


0 and a- = -Ba. 

5 4 


eliminating 4^2 equations (2.24) and (2.25), we have 
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the general solution of the equation (2.34) is 

^^4 ^ * -1/4 

= ^^9 - — G ^ 
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(2,35) 


where, a^ is the constant of integration. Equation (2.22) yields, 
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(2-36) 


The constant of integration a^ can be chosen to be zero 


and 


convenient values may be chosen for the arbitrary constants a^ , a^ 

and a^- . Thus, if we assume a. =1, a, = -1, = K and a^ = 0, 

o 1 4 o o y 

the expressions for 4^^ #2 ^ given by: 
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and 
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respectively. Finally, the governing equations reduce to 

f”’ + 3 ff" - 2f’^ - K f’ + d = 0 

o 




(2.39) 


Pr 


+ 3fd' = 0 


(2.40) 


with boundary conditions specified by the equations (2.11), (2.12) 

and (2.13). Thus, for steady free convection along a 

semi-infinite vertical wall with a constant surface temperature, 

similarity exists if K, the permeability parameter, varies as 
1/2 

X - > 

Case - II: Steady free convection with surface temperture 

varying as a power of a linear function of x. 

Since the flow is steady, a 2 = = a^ = 0. It is possible 

to obtain the general solution of equations (2.25) and (2.26) for 
the present case. 

A combination of (2.25) and (2.26) yields. 
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dx 


(2.41) 


where, « = The general solution of this equation is, 

o 4 

^2 = "10 

where, a^^^ is the constant of integration. Using (2.42), (2.25) 

and (2.26), the expressions for and obtained as. 




(.-1, 

Xi ajp 


= /(«-!) 


"lO ^ "ll a 


(«-l) X } 


10 


(2.43) 


Ue first discuss the case when « 1 . 

with the help of (2.43) and (2.44), 
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(2.44) 

is obtained from (2.22) 
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The constants a^^ , a^ , a^^^ and a^^^ are arbitrary. One can 
make a few simplifications here, without loss of generality. As 
can be seen from equation (2.45), a^^ can be chosen to be 1. By a 
similar reasoning, the value of a^ may be chosen leaving «, a^^^ 
and a., as arbitrary. Let us put n = - or « = From the 


definition 

of it is clear that 

= n- 1 , a 

- = n+3 and a- 

D / 

If 

from equations (2.25), (2.26) and 

(2.33) 

respect ively . 

The 

functions 4 ^ 
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form: 
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and the similarity equations reduce to 

f”’ + (n+3) ff” - 2(n+l) f’^ - K f’ + © = 0, (2.50) 

o 

i- d” + (n+3) {&’ - 4n f’© = 0, (2-51) 


which are solved numerically by the shooting method with boundary 
conditions (2.11), (2.12) and (2.13) for various values of n. 

The variation of physical variables within the boundary layer 
may now be expressed as: 
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It may be noted here that for n = = 0, the present case and 

case 1 are identical. It Is also observed from the present 
analysis that the numerical results are applicable for all values 
of a^j^ and n (positive, negative or zero), although all the cases 
do not correspond to physically realistic situations. 

The physical implication of non-zero a^ ^ is that the surface 
temperature at the leading edge of the plate is different from 
that in the ambient, a more realistic feature, since it is 
difficult physically to heat the plate without raising the 
temperature of its leading edge. For non-zero the velocity 

component u is not identically zero along x = 0 , but is zero only 
at y = 0 and y = oo. 

It can also be seen that when n = 1 and ®-4 becomes 

zero, indicating that the similarity variable, 17 is a function of 
y only and the plate surface temperature varies linearly with 
distance from the leading edge. The free stream velocity and the 
similarity variable for this situation are analogous to those of 
steady forced flow, in the neighbourhood of the stagnation 
point of a blunt-nosed cylinder . 

Yet another simplication arises when « = 1. The general 

solution for this case i« 

*2 ^10 ^ 1 ’ 

Thus, from equation (2.25) 
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The constants of inteeration a^^ and a^^ arbitrary. 

The similarity equations are, 


f”' + £f” - 2f’^ ~ ^ (2.58) 

d” + fe’ - 4f ’ d = 0 (2.59) 

rr 

with boundary conditions (2.11), (2.12) and (2.13). 

Unsteady Cases: 

For unsteady cases we find that the expression of 
permeability becomes a function of space as well as time. Ue omit 
the discussion of such situations as permeability is generally 
not a function of time in real physical situations. 

The above similarity equations derived for the generalized 
Darcy's model can be simplified to obtain the similarity form of 
equations for other formulations such as the Brinkman's model and 
the Darcy's model. A brief description of such simplifications is 
presented in the following sections. 

2^ ^ SIMILARITY SOLUTIONS FOR BRINKMAN MODEL 

Considering Brinkman model as the equation of motion, the 
generalized similarity equations (2.20) and (2.21) reduce to. 


^ - (a^n - f) « = 0. (2.61) 

For the steady case, a 2 = = 0 and a^, a^, a^ and a^ are, in 

general, non-zero constants. Using (2.26) and (2.28), one gets, 

G = A , 
w Z 


(2.62) 



The 


where, A ie the constant of inteeration and « = a /a . 
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expressions for ^ 2 * ^ ®ay now be obtained from the 

equations (2.22), (2.26), (2.28), (2.29) and (2.62) as follows: 

«-l 
4 ~4S 

^5^ ’ (2-63) 

3 

4 

^2 = (A^x + Ag) (2.64) 


3e 

= Aaj^ (A^x + Ag)^”*, ( 2 . 65 ) 

^ = ag Aj^(A^x + Ag)^"®" , (2-66) 

1/3 

where, A^ = (A/a^^) and A^ , A^ are the constants of integration. 

Let us first consider the case for € ?< 4 . On substituting n - 

.3« . . 

, one can write, 
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n-l 
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^^4* * ^5^ ’ 

n+3 

4>2 = CA^x + Ag)"~^ , 

= Aa^ (A^x + Ag)” , 

n-l 

^ = ag Aj (A^x + A_^)”^. 


(2.67) 

( 2 . 68 ) 

(2.69) 

(2.70) 


It is observed from 
transformation functions 
the same as in the case 
similarity equations are 
For Brinkman model, one 


equations (2.67) - (2.70), that the 

for steady flow using Brinkman model ard 
of the generalized Darcy's law. Thd 

, however, different for the two models, 
obtains: 


f”' - d - K f’ = 0, 

o 

d” + (3+n) f&' - 4n f'& = 0. 

For the case of « = 4, the expressions for , 4>2> 
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where, are the constants of integration. 

a. 5 SIMILARITY S(M-imC»IS FC»i DARCY* S LAW 

If the Darcy’s law is considered aa the iBooientuin equation the 
similarity equations (2.20) and (2.21) reduce to: 


f’ - K& & = 0, 


(2.77) 


^ d" - (a^Y) - a^ f) d’ - (a^ + a^ f’) & = 0, 


(2.78) 


where, a 
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For steady case a« = a. = 0. From the expressions (2.26) and 
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(2.28), one gets. 


^w = ^ ^2 ’ 


(2.80) 


where, A is the constant of integration and e = a^/a^. 

From equations (2.26), (2.28), (2.79) and (2.80), one obtains. 
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^w = * *12 ^*10* *11^' 
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Here K is taken as a 


constant and the above expressions are vali 
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^ 2. A.n, A., and A., are constants of integration. 
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On 


substituting n - , one gets, 

n-1 


^ ^12 ^^ 10 ^ 


(2.84) 
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“ (^lO* * »ll) 


= * *12 <*10* **ll’ 


n 


(2.85) 

( 2 . 86 ) 

I £ we choose a = 1, the above tranafornation functions are the 
same as those obtained by Johnson and Cheng [19] for a vertical 
plate embedded in a porous medium. The similarity equations 

(2.77) and (2.78) now reduce to the form: 

£• - kd = 0 (2.87) 

( 2 . 88 ) 


) fd' 

Pr 2 


n f’d = 0 


These equations have the same form as the similarity equatione 
obtained by Cheng and Minkowycz [13] and Johnson and Cheng [19). 

It may be noted, however, that the definition of the dimensionless 
temperature in the present work is different from those of [13, 
19 ], which gives rise to an explicit appearance of Pr and K in th< 

present formulation. 

The equations ( 2 . 81 )-(2 . 83 ) are valid for « ^ 2. For thf 

« = 2 one gets the exponential variation of the form, j 
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(2.91) 

are the constants of integration- 
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Sim & Itetrliuxl S^lut^lon 

The numerical aolutlon procedure for steady free convection wi 
vail temperature varying as a power function of x is presented 
her e . 


Let us consider the similarity equations (2.50) and (2.51) 
along with boundary conditions (2.11), (2.12) and (2.13). Ue 
employ the shooting method for solving these system of ordinary 
differential equations which reduces the solution of a boundary 
value problem to the iterative solution of an initial value 
problem. In this method the generalized similarity equations 
(2.50) and (2.51) are reformulated into a system of first order 
ordinary differential equations and solved. The solution 
procedure is briefly discussed below. Let, 
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The boundary conditions in terms of the new variables are: 



44 


At n = 0: 0, = 0, F^ = 1 , C2-97) 

= "^iHax** ^2 " ^ ■-■ ^2.98) 

The conditions at 17 = 0 are the prescribed initial conditions 


" ' «.a*= ’^2 ‘ O. F, ■= “■ 


The other required initial conditions (i.e. F- and F,. which ar 

J 3 


not specified at rj = 0 ) are assumed and the resulting initia 
value problem is solved by the Fourth-Order Runge Kutta method 
The assumed initial conditions are improved by the Newton-Raphso 
procedure until a solution which satisfies the known boundar 
conditions at far stream is obtained. The iterativ 
Newton-Raphson process can be set up as follows: 


Let, F-j rt = and F^ I „ = ol. and assume that 
3 It)= 0 1 .5«r)=0 2 


P-> 1 = (oi, , ot ) and F, I 

2 ' 17 = 1 ? 1 ^ 1 ’ 2 ' 4 ' 17=1 


= «2). 


max max 

The aim is now to solve for the roots of two numerically define 


functions Yj and Y^ which imply the satisfaction of the far strea 


conditions, as given below: 


S <^"l’ " 2 ^ 
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(2.99 

G 2 («^, « 2 ) 

- ^ 


= ^2 

(«!• " 7 ) 

= 0 . 

_ o._.. 

Using the generalized 

Newton- 

-Raphson scheme 

for finding the 

root 


of the non-linear algebraic system 

Y («) = 0, (2.101 
where, both the dependent variable Y and a are vectors. Knowin 
the value of « at the nth iteration, the (n+l)th iteration gues^ 
value can be obtained by the expression: 


Oi = a. - [j] -YCfit), 


where, [J] is the Jacobian ^ 


( 2.102 
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Since the functi 

are not known analytical 

the Jacobian the ^ v. j 

eueaa vai. , perturbed 

values of <» and are 

(for instance, as ot » ^ and 

“l - l.OOOl c, a„d »• = 1.0001 

nuiBerical di £ £ erem- 1 **“ ^ 

le „a»d to £l„d the 
a. 7 RESULTS AND Mscuag^ 


ta calculate 

aiy- ^ 


liable in the 


There are no detailed 

Q nuffierlcal results & 
literature for ^ in porous 

problems with free convecti 

media using eenerali^^. limiting 

Darcy's law. Therefore, * 

values of parameters K ^ ^ -,^a have been 

o *^d n, the present res^^ 

compared with those of „ond to free 

Datrach [26), which correeP 

convection over a vertical t i temperature. 

•»“ plate with constant wb-I^ 

Setting = 0 and n = 0 , the results 

tn the present £orffiulat ioo » 
for the situation studied w 

^ t>y Ostrach can be obtaii^^^'' 

The comparison betw^ [26] 

een our results and those ^ 

presented for n... .trslocity ar 

= 0.72. The 

_ 4 - 


been presi-***-^^ 

1-0 and Pr = 0.72. 

temperature distributi 

ons, are shown in tables 1 


ha< 


Table 1: Velocit-i, n, 

^l^y Distributi 


on for Pr = 1.0- 



0.2 

0.4 

1.0 

2.0 

3.0 

4.0 

5.0 

6.0 


0-10918 
0-18328 
0 . 25012 
0 . 13987 
D-05205 
D- 0 l 46 a 
0 - 00309 
0 - OQqoq 


0-1092 

0.1833 

0.2502 

0.1400 

0.0524 

0.0151 

0.0035 

0.0004 
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Table 2: 

Velocity 

Distribution 

for Pr = 0.72. 


Present Values 

Values in Ref. [26] 

0.2 

0.11581 



0.1159 

0.4 

0.19606 



0.1962 

1.0 

0.27543 



0.2759 

2.0 

0. 16S67 



0.1697 

3.0 

0.06443 



0.0661 

4.0 

0.01919 



0.0212 

5.0 

0.00415 



0.0060 

6.0 

0.00000 



0.0015 


Table 3: 

Temperature Distribution for Pr = 1.0. 

y) 

Present Values 

Values in Ref. [26] 

0.2 

0.88666 

0.8867 

0.4 

0.77421 

0.7742 

1.0 

0.46382 

0.4638 

2.0 

0.14220 

0.1422 

3.0 

0.03382 

0.0339 

4.0 

0.00712 

0.0072 

5.0 

0.00125 

0.0014 

6.0 

0.00000 

0.0002 


Table 4: 

Temperature 

Distribution for 

Pr = 

0.72. 

r> 

Present Values 

Values in 

Ref. 

[26] 

0.2 

0.89913 


0.8991 



0.4 

0.79886 


0.7989 



1.0 

0.51682 


0.5168 



2.0 

0.19435 


0.1945 



3.0 

0.05953 


0.0601 



4.0 

0.01618 


0.0170 



5.0 

0.00357 


0.0040 



6.0 

0.00000 


0.0009 




In figure 2, the velocity profile variation with 

permeability ia depicted. It ia aeen that with increasing 
(decreasing voidage) the velocity decreases in the boundary layer. 
This ia because of the increased resistance to flow at higher K^. 
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The velocity boundary layer is also thinner when the medium is 
less porous. For all values of , the velocity profile shows a 
typical maximum around y) = 1 as generally seen in natural 
convective boundary layers. 

In figure 3, the variation in temperature profile with is 
shown. It is seen that the temperature drop is leas steep when 
the voidage in the medium is small. This is because the rate of 
heat transfer is leas, when the flow is weak due to higher Darcy 
resistance. The thermal boundary layer thickness is only mildly 
influenced by the permeability parameter. Also, both the 
temperature and velocity profiles appear to be very sensitive to 
change in , at higher values of . Thus, at low permeabilities 
even a slight change could greatly affect the flow and heat 
transfer processes within the medium. 

The effect of Pr upon the velocity and temperature profiles 
are presented in figures 4 and S. In the purely Darcy model the 
effect of Prandtl number will not be felt because it does not make 
any explicit appearance in either energy or the momentum equation. 
It appears only in a combination of Rayleigh number for such a 
case. In the present work which focusses the non-Darcian effects, 
Prandtl number is seen to have strong influence upon both 
temperature and velocity profiles. 

The effect of n is shown in figures 6 and 7 on temperature 
and velocity distr ibut ions . 
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Figure 1 Physical geometry of the 
problem . 
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CHAPTER III 


CONJUGATE FREE CONVECTION FROM A VERTICAL FIN EMBEDDED 

IN A POROUS MEDIUM 

3.1 INTRC«>UCTION 

The problem of conjugate heat transfer from a downward 
projecting fin immersed in a saturated porous medium is important 
in applications such as the extraction of geothermal energy and 
the design of insulating systems for energy conservation. The 
problem discussed in the present chapter deals with a vertical 
plate fin situated in a saturated porous medium. The coupled heat 
transfer between the fin and the fluid flow is analysed in detail 
and a generalized non-Darcian flow model is used to describe 
natural convective flow within the porous medium. 

In the analysis of natural convective flows, it is possible 
to consider the fluid region alone exclusively, if the surrounding 
wall temperature is known. However, in several situations, 
especially for flows around slender bodies, a prescribed wall 
temperature distribution cannot be maintained. A majority of the 
available research studies dealing with heat transfer from 
submerged fins, assume that the convective heat transfer 
coefficient at the fin surface la uniform all along the fin. The 
fin heat conduction equation is then solved analytically using a 
uniform value of heat transfer coefficient [1,2]. However , the 
local fin heat transfer coefficient can, indeed, experience 
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eubatantlal variation 
nonuni form! ties in both 
fluid. 


aJong the fin 
the velocity and 


surface [3,4,.^] due tc 
teiBperature fields in th< 


Is necessary, therefore, to solve the conductive - 
convective heat tranaf^^ 

* a coupled problem and thereby 
usly analyse the flow and temperature fields on the 

the temperature distribution in the fin. The 

continuity of heat fiuv * 

and temperature at the solid-fluid 

interface ia demanded 

d, while solving the coupled heat transfer 

between the two phase* i 

aes- A problem of this kind is referred to as 

a conjugate haat 

Tha application ana thaory of heat tranafar by free 

convection fro. finned bodlaa haa bean atodlad for over half a 

century [6]. However & 

. a study of conjugate heat transfer for 

natural convective flou ti._ i 

was initiated by Lock and Gunn [3] in 

1968, who have studied rK., v j 

the boundary layer flow and heat transfer 

along a thin vertical fin rta w 

^^*1- They have shown that similarity 

solutions exist for certain types of surface temperature 
They also presented experimental results which agree 
theoretical predictions. Using a modified approach, Kuehn 
[7] also obtained similarity solutions for the problem, 
general solutions were obtained by Sparrow and Acharya fS] 
for natural convection about a vertical plate fin, with a portion 
the fin length heated by a plane heat source. Tolpadi and Kuehn 
[7] studied the problem of a vertical fin with a circular heat 
A numerical investigation of heat transfer between a 
rectangular fin and a combined forced and natural convective flow 
is presented by Sunden [ 9 ], 


Anderson and Be Jan [10,11] have 
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8tudi<^d the conjug,ate heat transfer problems for an impermeabl 
wall adjacent to a porous medium or a free fluid. However, thes 
solutions are applicable only to symmetrical problems in which th 
Oseen linearization can be used. 

The conjugate heat transfer for a fin embedded in a poroui 
medium has been investigated by Pop et al . [12]. They considerec 
the natural convective flow with a high Rayleigh number around a 
vertical fin. In another paper. Pop et al . [13] considered the 
mixed convection heat transfer along a vertical fin. They used 
similarity variables to transform the boundary layer equations and 
solved the resulting equations using a finite difference scheme 
based on central difference formulae. It is worthwhile to mention 
that all the above investigations treat the porous medium as 
purely Darcian. 

The purpose of the present work la to study the conjugate 
heat transfer problem along a fin embedded in a fluid saturated 
porous medium. A generalized non-Darclan model la invoked to 
describe the flow. Inertial and viscous terms are taken into 
account in the momentum equations. The governing equations are 
solved iteratively within the framework of boundary layer 
approximations by a highly-lmplicit finite difference scheme 
described by Hornbeck [14]. A variable mesh size is used in the 
direction normal to the fin surface, in order to obtain accurate 
flow and heat transfer characteristics inside the boundary layer. 

3. S GOVERNIMG EQLIATIONS 

A slender rectangular fin of thickness 2b and length L, as 
shown in figure 1, is considered. The fin is placed vertically in 
a fluid saturated porous medium and the base temperature of the 
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fin is taken to be prescribed. A laminar boundary layer flow Is 
assumed to exist around the fin surface. The fluid is 
incompr essibl e and the density is assumed to be constant except in 
the buoyancy force term according to the Boussinesq approximat ion. 
Viscous dissipation effects are neglected due to low flow 
velocities and one-dimensional heat conduction within the fin is 
considered since 2b << L. The momentum equation is described by 
the generalized Darcy ^ a law as proposed by Yamamoto and Iwamura. 
The governing equations for the fluid and the fin are therefore 
given by. 


35c ^ 

0, 


(3.1) 


Su 

= V + g^ (T-T ) - 

vti 

F" ^ 

o 

(3.2) 

u + 


" :;:7 ’ 

Oy 


(3.3) 

d T 

dx 

^f 

k b 
s 

r ^ ] =0. 

^ w 


(3.4) 

The boundary 

conditions are. 



y = 0: 

u = 

0, V = 0, T = T^, 



y eo: 

u 

0, T T 


(3..^) 

X = 0: 

u = 

dT 



X = L: 

T = 
w 

’■b- 



where, u 

and 

V are the horizontal 

and vertical 

veloci ty 


components, T is the fluid temperature and T^, and are the 

temperatures of the ambient fluid, the fin-surface and the 
fin-base respectively. Also, £ ia the acceleration due 

the kinematic viscosity, the coefficient of 


gravity, v, a are 
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volume expansion, and the thermal diffusivity of the fluid, is 
the permeability of the porous medium and , K are the thermal 
conduct ivitiee of the fluid and the fin. It may be noted that in 
the temperature boundary condition for the fluid at y = 0, the 
wall temperature is unknown and needs to be solved through the 
heat conduction equation for the fin. 

Introducing the following dimensionless quantities, 


T-T 

u = — , V = , & = 


b (Xi 


T -T 

^w T -T ’ * L ’ 
b ce> 

3 


(3.6) 


K g/S (T,.-T ) L" k,L 

y = y . K = , Gr = . CCP = ^ and Pr = - - 

^ L ^ s 


The governing equations (3.1) - (3.4) take the form: 


u 


dx 



Sv 

= 

0. 



a. tr 

Su 



^x 

T V 

^y 



^x 

+ V 

37 

= 

1 s^& 

¥r 2 

Sy 

& 

w 

.2 

+ CCP [ 


II 


+ Gr e - 


(3.7) 

(3.8) 

(3.9) 

(3-10) 


where, Gr is the Grashof number, Pr is the Prandtl number, K is 
the permeability parameter and CCP is the conduction~convection 
parameter. 

The equations (3.7) - (3.10) are to be solved with the| 

follo^jing boundary conditions, 

y=0: u=0,v=0,e= ©^(x), | 

y.>ao:u-^0,S.»0, 


X = 0 : u = 0 , d = 0 , 


d& 

dx 


w 


= 0 . 


(3.11) 


X = 1: d = 1. 

V 
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The boundary condition for d at x = 0 j s based on the 

V 

following argument. If the fin is very long and thin, the amount 
of heat which passes from the tip of the fin to the fluid la 
negligible, compared to the heat loss through the lateral Surface. 
In such a case the assumption of an adiabatic tip is justified and 
the temperature gradient in the fin at the tip is thus assumed to 
be zero [1,2]. 

It may be observed that In the foregoing description of the 
boundary conditions, no mention is made of the velocity condition 
at the base surface. This omission is required by the boundary 
layer equations. Once the conditions have been specified at the 
leading edge of the plate (x = 0), the solution yields the 

velocity field at all x > 0 and, therefore, does not permit the 
specification of velocity boundary condition at any x > 0. The 

non-accounting of hydrodynamic effects related to the presence of 
the base surface should not materially affect the fin heat 
transfer results provided that the fin length L is large compared 
with the boundary layer thickness at x = L [9]. 

3.3 METHOD OF SCAAJTICH 

Equations (3.7)-(3.10) are now expressed in finite difference 
form. A highly- implicit numerical scheme [14] is used in which 
the non-linear coefficients are taken at ( j + l)th level, where 
the subscripts j and k correspond to x and y directions, 
r eps ec t i vely . This form is necessitated by the zero velocity of 
the free stream, which will result in the classical implicit | 
formulation being inconsistent. The difference equations are | 


given as. 
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+ ~ ^3, k+l ^ '^j + l.k+l ~ ’^j + l,k _ „ 

Ax ■" = 0’ 


(3.12) 


u ^J^l,k ~ ^j,k ^ ^j^l,k^l ~ ^j^l,k-l 

5+1. k Ax j+l,k i(Ay) 


“j + l,k-i ■ ^'^j + l.k “j + l,k+l “j + l,k ^ 

2 K 

(Ay) 


®j + l,k' 
(3.13) 


^j+l.k ■ ®j.k ^ 

“j + l,k Ax '^j + l.k 


^j + l,k-i-l ^j + l,k-l 

i(Ay) 


d ~ 7A 4- A 

1 j + l.k-1 j + l,k j + l,k+l 

TF .. .2 

(^y) 


(3.14) 


Since the difference form is non-linear, an iterative scheme 
is required to obtain the solution. The difference equations 
given above can be written as, 

(i+1) (i+1) 


(i+l) V- 

“j + l,k+l ■ '‘j,k+l ^ ''j+l,k+l " ’^j+l,k _ „ 

_ - u, 

uCi+1) _ „ 

(i) “j + l,k '^j,k ^ „(i) “j + l,k+l j + l,k-l 

"j + l.k a 5E 5 + 1. k 2(Ay) 

j+l,k-l j+l,k “j+l,k+l 

= —-2 

(-iy) 

(i+i) 

5+1 -k ^ ^(i) 

— + Ur + 


u 


e e- 


^( 1 + 1 ) 


(i+l) ®j+l,k ^j.k ^ „(i+l) j+l.k+1 °j+l.k-l 

'^j + l.k Ax j+l,k 2(Ay) 


(3.15) 


(3.1<S) 


(i + l) _ 

1 ^j+l.k-1 J+l.k j+l,k+l 

TF 7~~1 


(Ay)‘ 


(3 


where, the auperscripts i and (i+l) indicate the values obtained 
at the ith and (i+l)th iterations, respectively. The Iterative^ 
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scheme is applied £or the above equations by solving one equation 
at a time in a ® manner. For instance, the momentum 
equation can be Handled in the following way. The simplified form 
of equation (3.16) £ 3 ^ 


(Ay 


(i) 

j + 1 .k 

2(Ay)- 
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(i +1 

j+l.k-1 
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(Ay)‘ 
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j+l,k ^ 1 

-A3r“ Y 





2(Ay) 
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(i + 1 ) 
j+ 1 ,k+l 
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(i + 1 ) 
j + 1 ,k 
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(i) 

^j,k 

Ax 


Gr & 


(i) 

j+ 1 ,k' 


(3.18) 


The matrix form of the above equation (after incorporating the 
velocity boundary conditions) is as follows. 
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1 0 
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j + 1. 2 
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A, B, c, 

0 s 
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D 
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(3-19) 
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where , 

B, 


(Ay)' 


rCi) 

■i + l,k 

^(Ay) 


u 


(i) 


(Ay)‘ 


(Ay)' 


* + i 

Ax 1( 


2 (Ay) 


(3-20) 


u 

. Gr , 

Ax J + 1 , k 


The discretized continuity equation (3.15) can be written as. 


j+l,k+l j+l,k 


Ay (i + 1) 

Ax ^'‘j + l,k+l “j,k+l^ * 


(3.21) 


Equation (3.21) is solved in an explicit fashion to obtain the 
velocity component v at all y locations starting from the 
wall. Now, the discretized energy equation (3.14) may be rearranged 
as , 


[~ 


Pr(Ay) 


(i + 1) 

""j-t-i.k e(i+i) + r 
2 2(Ay) ^ j+l,k-l '■ 
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Pr(Ay)' 


(i + 1) 
j-^l>k 
Ax 


^ 1+1. k 


+ [- 


(i + 1) 

. ^j+l.k , ^(i+1) 

2 2(Ay) ^ j+l,k+l 


j+l,k j,k 
Ax 


Pr(Ay) 

The above equation can be ptit in the matrix form. 


( 3 . 22 ) 


after 


incorporating the temperature boundary conditions as follows. 
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(3.23) 


(3.24) 


Thft fin heat conduction equation (3.10) nay be written in the 
difference form as. 


. 

J-1 J J + l 

(Ax)^’ 


= - CCP f — 1 

^ I Vj,0) 


(3-25) 


where d'f is the fin-surface temperature at location j. The above 
J 

equations can be rearranged in the more convenient form. 



bY 


j_l - 2^'? ^'j + 1 = ■ f I 7 ^(j,0). 


(3.26) 


After incorporating the boundary conditiona for the fin, the 
matrix equation corresponding to equation (3.26) can be obtained 

easily. 

Due to the non-linearity and the mutual coupling between the 
equations (3.19), (3.21), (3.23) and (3.26), an Iterative 

solution needs to be developed. The iterative process is complex 
since it involves an overall iteation loop between the fluid and 
the fin equations and a sub-iteration loop for the fluid phase 
equations. The basic steps of the procedure are summarised below: 

(a) An initial guess for the fin temperature S^(x) is estimated. 

2 

A parabolic variation, say ^^(x) = x was found to be a good 
initial guess. 

(b) The flow field and convective heat transfer equations in the 
fluid are solved with the guessed fin temperature 
distribution as the wall boundary condition. 

(c) From the calculated fluid temperature field the heat flux 

is estimated at every x-position. 

'-^y-'w 

(d) The fin heat conduction equation is solved with the 

prescription of the value of (||)y as obtained above. The 

for the fin temperature distribution ^„(x) is 
tl0W ^U0SS I w 

evaluated and steps (b-d) are repeated until convergence is 
btalned Uhen acceptable convergence of the above mentioned 
iterative process is achieved and the conditions of 
continuity of heat flux and temperature at the fluid-solid 

interface are 


satisfied, all the relevant heat transfer 


characteristics can 


then be calculated. 
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The equations have been solved by using a constant step-size 
in X direction and a variable mesh size in the y-direction. A 
finer mesh size is used in the region of rapid variation i.e. near 
the fin and a relatively coarse mesh is used for the region of 
slower variation. 

There are many advantages of using the variable mesh 
technique over that of a uniform mesh-size over the entire region. 
Firstly, with a variable mesh-size, the number of simultaneous 
equations to be solved can be reduced, which is very important for 
an implicit formulation. Since the time required to solve the 
simultaneous linear equations increases with the number of 
equations, a considerable saving in computational time can be 
achieved. In addition, the round-off error accumulated while 
solving the large number of simultaneous equations is held to a 
minimum. Furthermore, since the boundary layer has zero thickness 
at the tip of the fin, a uniform discretization will have very few 
grid points near the fin tip. Hence the solution will not be 
accurate near the tip and these inaccuracies will be transmitted 
along the full length of the fin. In the present work a 
three-stage mesh size variation has been considered and due care 
has been taken to provide sufficient number of grid points within 
the boundary layer region near the fin tip. 

3.-4 RESULTS AND DISCUSSION 

The temperature distribution in the fin and the porous medium 
are obtained along with the velocity distribution in the fluid, 
for a wide range of governing parameters of the problem, namely 
Gr, CCP, Pr and K. A constant step size of Ax = 0.002 is used in 
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the x-directlon while a three-stage variable mesh size is used in 
the y-direction. The step size Ay varies as follows. 

Ay - 0.002 for the first 51 steps from the wall. 

Ay = 0.005 for the next 130 steps and 

Ay = 0.01 for the last 6 steps. 

The total y distance covered by these steps, (y ), approximates 
the thickness of the boundary layer in the present calculations. 

To establish the correctness of the present work, the results 
obtained by the finite difference method are compared with 
analytical results for a limiting case. The situation considered 
for this purpose corresponds to a low Grashof number for a free 
fluid (K oo) . For small Gr, the temperature solution in the fluid 
around the fin is dominated by one-dimensional conduction normal 
to the fin and the horizontal and vertical components of fluid 
velocity are very small so that the convective effects can be 
neglected from the energy equation. Uith these assumptions the 
energy equation reduces to. 



The general solution of the above 
given boundary conditions is 


(3.27) 

equation which satisfies the 


& = d^(x) (1 


max 


), 


where, & (x) = and m^ - CCP/Y , for CCP * 0. 

’ w cosh m »ax 

= 0, there is no heat loss from the fin and the entire 
taken to be at base temperature. So, one can set (x) 


(3.28) 

For CCP 
fin can 
= 1 for 


CCP = 0. 


Neglecting the inertial and the Darcy resistance terms 
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equation of motion can be written in a simplified form as 

^-4 + Gr d = 0. 
dy^ 


(3.29) 


The solution of equation (3.29) with the given velocity boundary 
conditions is, 


u = Gr 



(3.30) 


In Table 1 the results corresponding to the above described 
limiting case, obtained by finite difference method are tabulated 
with the exact analytical values. 

Table 1 

3 

Temperature and velocity distribution for K ot> (10 ), CCP = 0, 

Gr = 1.0, y =0.81. 

’ 'max 


Temperature Horizontal velocity 



Nunierical 

results 

Analytical 

results 

Numerical 

results 

Analytical 

results 

0.0 

1.0000 

1.0000 

0.0000 

0.0000 

0.1 

0.8780 

0.8765 

0.0225 

0.0222 

0.2 

0.7560 

0.7531 

0.0363 

0.0356 

0.3 

0.6340 

0.6296 

0.0425 

0.0416 

0.4 

0.5120 

0.5062 

0.0423 

0.0412 

0.5 

0.3901 

0.3827 

0.0371 

0.0357 

0.6 

0.2681 

0.2593 

0.0279 

0.0264 

0.7 

0.1462 

0.1358 

0.0160 

0.0146 


The j^ood agreeiaeiit between 

the numerical 

and the exact 

analytical solution 

indicates the 

correctness 

o£ the results 


obtained by the finite difference method. 
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The wall temperature variation with x ia plotted for varioua 
Gr in figure 2. The analytical result obtained for small Gr is 
shown by a dotted line. It is found that as Gr decreases the 
temperature profile approaches the hyperbolic variation as 
suggested by the analytical solution (CCP * 0 and Gr low). As Gr 
increases more heat is removed from the fin by fluid convection as 
indicated by the slope of wall temperature. In figure 3, the wall 
temperature variation for various permeability parameters is 
shown. It is seen that as the permeability parameter K approaches 
the value of zero, the wall temperature distribution moves closer 
to the analytical solution for pure conduction in fluid. This is 
to be expected, since for small K, the Darcy resistance to flow ia 
very high which reduces the velocity of the fluid. The heat 
convection effect, therefore, becomes very small even at large Gr 
as K approaches zero. thus for increasing K the fin is cooled 
more. 

In figure 4, the effect of the conduction-convection 
parameter, CCP, upon the wall temperature distribution la 
depicted. For each value of CCP, the numerical solution and the 
corresponding pure conduction solution have been plotted. From the 
figure, it is observed that wall temperature increases with CCP. 
This trend is easily understandable since higher CCP leads to more 
coupling between the solid and fluid phases and hence more heat 
loss from the fin. Both the analytical and numerical aolutioni 
exhibit similar variation with CCP, although a close match 
between the numerical and its corresponding analytical solution is 
not to be expected on account of the high Gr value used for the 


numerical solutions. 



12 


Ttue variation of wall temperature diatribution vLthi Pc', 
Pc'amdtl number ie shown in figure 5. It ia seen that the fin 

cooHi^ft is more effective at higher Prandtl number due to a 
thiitnn er thermal boundary la yer . 

1ft figure 6 the fl^id temperature variation in jc-d ir ec tl on. Is 
ahiowft f er various Graahof nmanbers at a fixed value of y. The 
fh uJI-d temperature has a sloiHer variation as the fin temperature 
vL tha rreapect to x. As Graahof number increases the temperatu re of 
the fmuJd decreases. This trend is similar to that of the fin 

tesippec'aitu.re. Although the total amount of heat removed from the 

fl n ies larger at higher Graah-of number, the increase In the heat 
-tran-sEetr coefficient is more rapid resulting Lji a smaller fluid 
temp er-afeure. Figure 7 and d show the variation of the velocity 
ti el d and the temperature field in the y-direction, at ■various 
>r- lo catLoms . As expected. Cor fixed values of other parametera , 
the -valtmes of velocity and temperature increase with x. 'The 

i-ftcrr eases In the velocity values with x indicates that natuaral 

ccociv ec tLon becomes more and aore vigorous as fluid flows adjacent 

to t»ie fin in x-direct ion. This is a direct consequence of the 

Lncreaaetd heat penetration Into the fluid at higher valtiea of x_ 

Ian figures 9, 10, tl and 12, the local Husselt noalner 

vatriat io n with the governing parameters Gr , K, CCP and Pr are 
tthovtn. The Nusselt number increases with Gr , K and Pr and it 
deccreaaiea slightly with CCB. Since the heat loss from tha flrx if* 
related to the slope of the wall temperature variatLora in 

jc-ditreti^n and the Nuaeelt truaber characteriaea the heat loas ftoaa 
tha £ic\, the observed variation of Nusselt number with vrar'ious 
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parameters can be explained by considering the slope of in 
figures 2-5. 

In figure 13, the wall temperature profile predicted by the 
present analysis is compared with that of Pop et al . [13] who have 
considered Darcy's model for flow. For smaller Grashof number the 
difference is not that significant. However, as expected, the 
non-Darcy effects play an important role at higher Grashof 


numbers . 
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Figure 1 Fin Geometry 



Figure 2 Wall Temperature versus x . 




Figure 3 Wall Temperature versus x . 




' I 1 1 I I 1 I 1 1 1 1 

0 0.5 1.0 


X ► 


Figure 4 Wall Temperature versus x . 



Figure 5 Wall Temperature versus x . 
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Figure 7 Velocity distribution versus y. 
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Figure 9 Nussclt Number versus x . 
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Figure 11 Nusselt Number versus x. 
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Figure 12 Nussclt Number versus x . 




Figure 13 Comparison with [13]. 



CHAPTER IV 


UNSTEADY NATURAL CONVECTION IN AN ENCLOSURE 

4. 1 I NTRCMMJCTI ON 

Natural convective heat transfer in enclosed porous media is 
a fundamental problem that pertains to many practical engineering 
situations. Applications include thermal Insulation systems, 
nuclear waste storage devices and porous heat pipes, just to name 
a few. 

The present study deals with the problem of natural 
convection in an enclosure which surrounds a saturated porous 
medium. The enclosure is made up of two isothermal vertical 
surfaces and two adiabatic horizontal surfaces. The walls of the 
enclosure are impermeable. The need for studying such problems 
had arisen from the common application such as the utilization of 
air gaps for effective insulation in the construction of 
buildings. It is a common practice to build brick walla separated 
by an unventilated air gap of few inches while constructing 
dwellings. For many years, heating engineers have been concerned 
with the calculation of heat transfer across this air gap. The 
problem of this type was first formulated in two-dimensional form 
by Batchelor in 1954 [1]. 

An appreciable insulating effect may be achieved by placing a 
porous material (fibre glass) in the unventilated air gap between 
vertical walls. This is due to the fact that multicellular 
convection does not occur in this case as shown theoretically by 
Gill [2] for a porous slab of finite height. Experimental 



observation, .ad. by Klar.f.ld [3) and Borl.. Co.barnou. [4], 

and the numerical predictiona of Chan et al rti 

L^J, ^180 ifldiCAllft 

that a unicellular motion exiata inaide the poroua gap. 

The problem of thermal convection in a confined poroua medium 
has also been of some interest to geophyaiclsta . In many 

naturally occurring eituatlons, fluid motion results from the 
variation in the buoyancy force, which is a consequence of 
temperature gradient in the fluid. Thia fact has led to the 
theoretical and experimental studies of natural convective flow 
inside enclosures of various shapes filled with saturated poroua 
media. A substantial amount of work in this regard has been done 
on the rectangular enclosure, particularly when either the 
vertical or the horizontal walla are maintained at constant 
temperature, while the other walla are insulated. Convective heat 
transfer in a rectangular poroua enclosure, whose vertical walls 
are maintained at two different temperatures and horizontal walls 
are insulated, is a problem of fundamental interest in this area. 
This problem has been studied numerically by Chan et al. [5], 
Vlasuk [6], Holst and Aziz [7], Rankvall [8], Burns et al. [9], 
Hlckox and Gartling [10], Prasad and Kulackl [11,12]. The 
boundary layer analysis of the problem has been presented by 
Ueber [13], Ualker and Homsy [14], and Bejan [15], while Simpkins 
and Blythe [16] have presented integral solutions. Approximate 
solutions for the same problem have been obtained by Ualker and 
Homsy [14], and Rejan and Tien [17]. Important experimental 
results are reported by Klarsfeld [3], Bories and Combarnous [4] 
and Seki et al. [18]. 



Moat of the above analytical and nuaerlcal 




Darcy's law as the laomentuffi equation and only a few researchers 
have considered non-Darcian flow behaviour Jn porous media. Chan 
et al. 15] have adopted the Brinkman model which Includes viscous 
diffusion terms in the fluid. They Inferred that the inclusion of 
viscous force term and no-alip condition on the Impermeable 
enclosure valla give virtually the same result as a pure Darcy 
analysis for any reasonable value of Darcy number. Tong and 
Subramanlan {19] gave a boundary layer analysis for natural 
convection in a vertical porous enclosure on the basis of Brinkman 
model. A numerical solution for a high value of the modified 
Rayleigh number was given by Tong et al . [20] using the Brinkman 
model. Inertial effects on buoyancy driven flow and heat transfer 
in a vertical porous cavity were examined by Prasad et al . [21] 
using the Forchhelmer-extension of Darcy equation. 

It is worthwhile to mention that generally, researchers have 
concentrated mainly on steady state natural convective flow in 
rectangular cavities with both vertical walls at uniform but 
different temperatures. Comprehensive reviews of the existing 
literature has been presented by Ostrach [22, 23] and Catton [24] 
for the free fluid case. However, transient development of heat 
transfer and flow may be of importance in many natural convective 
problems [25] in the porous media. 

The objective of the work presented here is to study the 
two-dimensional unsteady natural convection inside a rectangular 
porous enclosure saturated with an incompressible fluid. The 
enclosure consists of two isothermal vertical fcj^lls and two 
adiabatic horizontal walls. The flow and temperature field 


eolutiona are obtained using the finite difference technique. The 
upwind differencing procedure is applied for convection teras 
along with central differences for diffusion terms, for spatially 
discretizing the governing equations. The timewise discretization 
has been performed using a fully implicit scheme. The resulting 
equations are solved using the line-by-line technique coupled with 
the tridiagonal matrix-solver algorithm. 

A parametric study is presented for various governing 
parameters of the problem such as the Rayleigh number, Ra, the 

Darcy number. Da, and the Prandtl number, Pr. Streamlines and 

isotherms are plotted for different times, to analyse the flow 
patterns inside the enclosure. The transient variation of the 
average Nusselt number at the hot wall has also been obtained in 

order to analyse the rate of heat transfer between the wall and 

the fluid trapped inside the porous matrix. 

4.2; 60VERN11«S EOUATIONS 

A rectangular, two-dimensional enclosure with isothermal 
vertical and adiabatic horizontal walls is considered (figure 1) 
1. Initially the fluid is stationary and has the same temperature 
T as that of the walls. At time 1 = 0, the wall at x = 0 is 
suddenly brought to a temperature (T^ > maintained 

at this temperature for t, > 0. Thus, the heat conduction and 

thermal convection processes are automatically set-up in the 
fluid. The following assumptions are invoked to simplify the 
problem formulation: 

a) The Boussinesq approximation is valid, so that the variation 

in density can be neglected except in the buoyancy term of 

This is true if the temperature 


the momentum equation. 



difference between the iaothermal walla and the height of the 
cavity are not very large. 

b) The flow is two-dimenaional and laminar. 

c) The fluid under conaideration ia Newtonian and all propertiea 
except density are conatanta. 

d) The compression work and viacoua dissipation are negligible. 

e) The boundary walls of the enclosure are impermeable. 

f) The enclosure is filled with a porous material which ia 
saturated with fluid. 

Thus the mass, momentum and energy equations for the flow 
inside the porous medium are. 
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The initial and boundary conditions of the problem are: 
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= B for 0 ^ X ^ A. 


where, u and v denote the velocity components in x and y 

directions respectively, P is the dynamic pressure, T 

fluid tumpurature. K la tha per«a.billty of th. ..dl«« «.d « i« 



the acceleration due to sravity. The properties p, « and f) 

denote viscosity, thermal dlffusivity and coefficient of thermal 

expansion respectively and p is the density. Also, T and T 

’ H c 

denote the temperatures of hot and cold walls and the constants A 
and B represent the horizontal and vertical dimensions of the 
enclosure . 

The following dimensionless quantities are Introduced, 
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where, Ra and Da are the Rayleigh number and the Darcy number 
respectively. As is the aspect ratio and Pr is the Prandtl 
number . 

The equations (4.1) - (4.4) may be written in their 

dimensionless form as, 

(4-7) 
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(4.10) 


The equations (4.7) - (4.10) describing thermal convective 
flow inside the porous enclosure, are a set of simultaneous 

non-linear and elliptic partial differential equations. Complete 

analytical solutions for these equations are not possible and one 



needs to take recourse to numerical solution techniques. In the 
present study the finite difference method based on stream 
function - vorticity flow formulation has been adopted for 
convenience in numerical simulation. 


Eliminating p from equations (4.8) and (4.9), and introducing 
stream function y* and vorticity < 0 , as 


u 







and 


(4.11) 

(4.12) 


the following set of equations is obtained: 
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(4.14) 

(4.15) 


The initial and boundary conditions to be satisfied by y>, o and B 
are given as follows: 
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(4.16) 

0, V " 0, ^ . 

0, 


0.^.0,^ = 

0, 



T > 0. r, - 0= 3^ 

T > 0. „ = As: ^ = 0, w = 0, ^ = 0. 

The difficultly aaetociated with the deteriainat ion of preaaure field 
has led to the stream function - vorticity method, which 
eliminates pressure from the governing equations. Moreover, in 
this formulation, instead of dealing with continuity equation and 



two mooientua equations. 


one needs to solve only two second 


order 


equations . 

4.3 METHOD OF SOLUTIOM 

The governing equations (4.13) - (4.15) for unsteady free 

convection are coupled, non-linear partial differential equations, 
the solution of which eventually converges to the steady state 
solution of the problem. A marching procedure in time is employed 
to determine the velocity and temperature distributions at various 
times. Using an implicit scheme the governing equations may be 
written in the difference form as follows: 
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where, i and j are subscripts corresponding to ( and yj directions, 
41.0. A,.. = - V-?].! 

^nd the superscript k represents the tiiae level. 



To obtain the difference equations of (^-17) and (4.18), 


Upwind-differences for the convective terms and central 
differences for other terns have been employed. For obtaining 
equation (4.19), only central differencing has been used. In the 
Upwind differencing scheme one uses the backward difference 
expression for a derivative term when velocity 1*^ concerned 

direction is positive and a forward difference expression when it 
is negative. It is evident from above difference equations that 
they satisfy the "positive coefficient rule" required for 
stability [261 and therefore, the solution from the present 

*iunerical formulation is stable for all step-sizes and Ayj 

_ , Tjy.ii'fliffb em e has be en 

Spatial discretization. Moreover, since implicit a 

n - j /- . 1.1 ani 1 nnA are stable for all 


^Sed for timewise marching, the solutions are 
step sizes also [27]. The truncation error in 


the final solution 


0 (A^: ,Ar?, At) . 

«ived by line-by-line 

The resulting difference equations are s 

1 auiti-dinensional 

Method which is a popular technique for solving 
r. A solution 
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^^erations, for each time step. During the horiz 
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aoat recently available valuea. boundary 


^ .v from 

^<ial equations for (1, j) snd (H+ , £ or each J 

puf set up 

^tiditions, a pseudo tridiagonal system wi 
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which can be solved by tridiagonal aleoriths. After sweeping all 
the horizontal lines except those with prescribed boundary 
conditions, if any, one can perform vertical sweeps in a similar 
fashion by transferring terms corresponding to (i + l,j) and (i-1 j) 
to right-hand side for each i. The line-by-line sweep can be done 
for the horizontal and vertical directions alternately until 
convergence for all grid points occurs. 

The equations (4.17) - (4.19) can be rewritten in the 

triadiagonal form as follows, for applying the line-by-line 
procedure. Considering equation (4.17), at a point (i, j). 

For Horizontal Swoop : 

Case 1 : > 0 and Ay»^ > 0 

k+ 1 

(-(Ay^l - Pr) J + * 4rjPr + 4r2pr + r^ 

♦ 2 An Pr . (4.20a) 

= 2rj Pr * (|Ay.| a Zr^Pr) * rj 


-An RaPr - «^!i_ 5 ). 
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- An Ra Pr (d.^^ . - j). 



98 


Caae 3 : > 0 and Av'^ < 0 


C-jAv'^j 2 r 2 Pr) ^ + (jAv<.| + jAy/jl + 4rjPr + Ar^Pr 

(4.20c) 
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Case 4 : Ay^ < 0 and Ay^^ > 0 
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For Vertical S«reep 

Case 1 : A^'^ > 0 and Ay^ > 0 
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At> Ra Pr 
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Case 3_£_ Av'^ > 0 and < 0 
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The boundary conditions for are obtained as follows: 

^2 ^2 

From the definition of vorticity, oi = — y + — y. Using 


the 


no-slip conditions on the boundaries in the above equations for u 
one can obtain the following expressions 


„ = £!| at ? . 0, 
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O w ^ 
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itK »r, 


(4.21) 


The surface vorticity values can now be evaluated through the 
finite differences form of equation (4.21) for the concerned 


boundary . 


TOO 


(i) At = 0: 


Applying central difference acheme to one obtains 
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(4.22c) 
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Let us now consider the line-by-line solution of eq. (4.18). 

For Horizontal Swoop 

Case 1: Av . > 0 and Av. > 0 
J I 

(-lAyi^l - Zr^) + lAy-jj + 4rj + 4r2 + r^) 

- 2r 
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.k + 1 


k+1 
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Case 2 : Ay^ . < 0 and . < 0 
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Caaft 3 = > 0 and < 0 


,k+l 


(-|Av.| - 2r^) KV'-I + 4rj + Av^ + r,) ©?- 


k+1 

3" ^i.j 


« x.k +1 

2 r* & 

^^2 i+l.j 


(4.23c) 


.k+1 


(|A»..| - 2rj) ♦ 2rj ^ . , 


Caae 4 : Ay'^ < 0 and Ay^. > 0 


-2r„ 


k+1 


r-j.j * * i-^ji * ‘■'i * '■^2 ^ ■'3> *tj * 


k + 1 


C-|AVjl - 2 - 2 ) 

= 2r, * Cliv-jl * 2r3) ♦ ^3 j. 

For Vortical Swoop 


(4.23d) 


Case 1: Ay. > 0 and Ay. > 0 
J ^ 


o x.lc+1 

-2r. d. , , + 


1 -i,i-l * “'1 * “^2 * -3- >,j 


r,) * 


(-|A»..| -2r3) 


(4.23e) 


.k+1 


.k+1 


k+1 


ClAVjh 2r3) * 2r3 * Tj a . _ , . 

Case 2: Ay. < 0 and Ay. < 0 

— j j. 

(-lAVjl - 2 r 3 ) * lAv'jl + *>-3 ♦ 4r3 * .Tj) sr: 

- 2r, a*?*' , (*-“0 

1 1 , j + 1 


Caae 3 : Ay^ > 0 and Ay^ < 0 


k+1 


(-|Ay^| - 2r^) + ( lev'll * ^^2 * ®i,j 


- 2r, &. 

1 i , j+1 


(4.23e) 
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Caae ^ Av'j < 0 and > 0 


-2r, e 


k+1 


1 ®i,j-l * 'l^'"il * l^*'jl ♦ + 4rj, ♦ Tj) 8^* 


k+1 
j 




.k+1 


(4.23h) 


= 2*^2 ^i-l.j ^ Cl^V'.|+2r2) . + r. d 


k 

3 “i,j 


The line-by-line procedure ia applied to equation (4.19) aa 
followa : 

For Horizontal S^woep 

= -r5-i - <5.1 * <5- 

For Vortical Swoop 


2 k+1 2^ k+1 , 2 k+1 

r V'- .,-2(l+r)v'. .+r yi. ... 

1 , J-1 -1,1 1, J+1 

k+1 k+1 ^ k+1 .. 

' - n-t.j - »'ui,l * “l.j ■ 

where, r = ( ^ )^. 

A auitab] e under-relaxation factor ia used for each equation to 
increaae the stability of the iteration acheoe. 

An important physical characteristic of the problem la the 
rate of heat transfer from the hot wall to the fluid. In the 

present problem, the Nuaselt number on the hot wall of the 
enclosure ia calculated by 

„„ = ^1 (4.25) 


— 


y )=0 
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Th« wall gradient i. calculated u.l„g . three-point for.ule. 
After obtaining the gradlenta at each grid point lying on the hot 
wall, the average Nu.aelt nuab.r 1. calculated by Slap.on'. 
integration fornula. 

The computational procedure adopted for obtaining the overall 
solution of the problem may be eummarised as follows, 

(a) Initially at t = 0 . set d = 0 and w = 0 , as given by the 
initial conditions of the problem. 

(b) Set T = T + At , where At is the time step. 

(c) Guess the values of 6 and y at the new time step. 

(d) Get the nodal values of B from equation (4,18). 

(e) Calculate ^ from equation (4.17), using the d-values obtained 
in step (d) . 

(f) Obtain the values of by substituting the values of <•> from 
step (e), in equation (4.19). 

(g) Iterate steps (d - f) for a few times. 

(h) Uith the current values of B, ui and y, proceed to the next 
time step (begin at Step (b)). 

(i) Repeat the above process is repeated until steady state is 
reached. 

(j) Calculate the local Nusselt number using a three-point 
formula for the temperature gradient at the hot wall and then 
evaluate the average Nusselt number by Simpon’s 1/3 rule. 

A. A RESULTS AND DISCUSS! CM 

Results have been presented for the temperature field, the 
flow pattern and the average Nusselt number at the hot wall, at 
different times and for various parameters. 
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ord.r to v.rify th. .ocutacy, o, th. n„U. dUfaoaoc. 

dl.crotlzatloo In apaca, th. t.a.u, ealcul.t.d oal„e 

different atep-alnea for the encloaur. aapect ratio of 10 . the 
value of the ateady atate Nuaaelt nuaher „aa found to have „o 
aignif leant changa for the meshes of 60 x 60, 80 x 80 and 100 x 
100. Thus, it may be concluded that the finite difference spatial 
discretization is adequate. In the present analysis the 80 x 80 
grid system has been used and a fixed aspect ratio of 10 has been 
considered, for all the calculations. 

The flow patterns inside the enclosure at small and large 

times after the start of thermal convection inside the enclosure 

are shown in figures (2-5) for various sets of parameters. It is 

clearly seen from the plots that unicellular convection occurs for 

the chosen aapect ratio and range of parameters. Similar results 

have been reported for a vertical layer by iClarsfeld [3]. For 

smaller time the circulation vortex is closer to the hot wall 

since the convection starts at time t = 0, near the wall. Also the 

close spacing of streamlines adjacent to the hot wall for small 

times indicate that a boundary layer exists there. For later times 

due to vorticity diffusion into the interior portion of the 

enclosure, a circulation vortex which is almost symmetric about 

the midplane develops. The total dimensionless time to approach 

steady state depends on the modified Rayleigh number, Ra^ (= Ra x 

Da). It is observed that ateady state is attained at a smaller 

value of time for small Ra . 

o 

At a fixed Darcy number the boundary layer near hot wall 
becomes steeper as Rayleigh number is Increased. This can 
easily understood from the fact that if the temperature poten 



or fluid vl.oo.lt, d.or...... th. bound.,, l.,.o „ni 

b.co.. .t.ep.r. It 1. . 1,0 ...„ 

function. In Fleur,. (2-5) th.t th. .tr.neth of th. ciroul.tlon 
vort.x Incr..... with R.yl.leh nuBb.r, For th. r«>g. „f R.yl.igh 
numbers considered, this increase in vortex strength is almost 
propotional to the Rayleigh number. The shapes of circulation 
vortices indicate that near the corners of the enclosure the flow 
is virtually stagnant. For a fixed Rayleigh number, when Darcy 
number is reduced the strength of the vortex decreases and the 
boundary layer becomes less steep due to high resistance of the 
porous matrix. From the above trends it maybe concluded that for 
thermal convection problems in a porous medium, the product Ra x 
Da is more relevant than the Rayleigh number or Darcy number, 
individually. This feature has also been observed by Chan et al. 
[5] earlier. 

The shapes of isotherms at small and large times are shown in 
Figures (6-9) for various values of the parameters. It is 
observed that for initial times the isotherms are closer to the 
hot wall since heat transfer la affected by contact with this wall 
at starting time, t * 0. For larger times the isotherms are more 
evenly distributed due to the heat transfer mechanisms of 
conduction and thermal convection. The shapes of the isotherms are 
almost vertical straight lines for small times, indicating that 
the heat transfer within the enclosure is dominated by conduction 
in the beginning, before thermal convection gains strength. In 
general, the isotherms are curved towards the hot wall at the 
bottom of the enclosure and curved away from the hot wall at top 
in exact correspondence with the direction of flow within 



enclosure. The effect of modified Ravl«.iflK 

Kayieigh number upon the 

isotherm patterns is as follows. At low values of Ra , the steady 
state isotherm pattern also consists of evenly-spaced, 
nearly-straight lines, indicating that the effects of convection 
are minimal and the heat transfer is dominated by conduction. For 
larger values of , the Isotherm patterns become more and more 
distorted at large times, when thermal convection is fully 
set-up. This is clearly shown by the shift of the isotherms 
towards the hot wall at the bottom of the enclosure and away from 
the hot wall at the top. Thus, the shape of the isotherms is also 
determined by the modified Rayleigh number rather than Ra and Da 
individually, as seen in the case of flow patterns earlier. 

In Figures (10, 11) the transient variation of the average 
Nusselt number at the hot wall is plotted for different Rayleigh 
and Darcy numbers. The effect of increasing Rayleigh number or 
increasing Darcy number is similar upon the Nusselt number 
variation, confirming the fact that Nusselt number is also 
determined by Ra x Da. At time x = 0, the Nusselt number is 
infinitely large due to an infinite gradient existing when the hot 
wall is suddenly brought Into contact with the cold fluid. The 
Nusselt number decreases gradually with time (at low Ra or Da) 
upto the steady state value due to the continuous penetration of 
heat into the porous medium and the consequent decrease in the 
value of the wall temperature gradient. At high Ra or Da, the 
Nusselt number increases briefly after the Initial sharp decrease 
when the wall comes into contact with the fluid. This increase in 
the Nusselt number may be attributed to the existence of strong 
convection currents adjacent to the hot wall when heat has 



IKJI 


penetrated completely into the porous laver tvs^ ^ 

xayer. ifte Strong buoyancy 

effects adjacent to the hot wall at high Ra (aided by the low 
resistance from the solid matrix at high Da) will generate a 
vortex flow of high strength adjacent to the hot wall for small 
times. This strong convection, in turn, tends to increase the 
wall temperature gradient and consequently an increase in the 
Nusselt number at the hot wall. Uhen steady state la achieved, 
heat penetration through the fluid is complete in the entire 
enclosure and the strength of convection existing adjacent to the 
wall also slowly decreases due to the fluid mixing. This leads to 
the decrease in the wall temperature gradient and in the Nusselt 
number value. 

A comparison between the present results and those of Bejan 

[15] are presented in Figure 12. The graphs have been plotted for 

Nusselt number versus modified Rayleigh number for the aspect 

ratio of 10. Ueber’s [13] results are also shown in the graph. 

It is observed that the present results are very close to that of 

Bejan and Ueber for lower values of Rayleigh number. However, for 

larger values of Ra the results are a bit deviated from Bejan’s 

o 

results which are for the Darcy’s law. This difference may be 
attributed to the inertial and viscous effects taken into account 
in the generalized Darcy’s model of the present study. 
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(i) Temperature boundary (ii) Finite difference 

conditions for the problem mesh 


Figure 1 Geometry of the system under 
consideration. 
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CHAPTER V 


FREE CONVECTION FLOW ALONG A MOVING VERTICAL WALL IN A POROUS 

MEDIUM 

B. 1 IhITfiODUCTJON 

Heat transfer from a vertical surface embedded in a saturated 
porous medium finds application in many branches of technology 
such as the power production from geothermal energy and the 
design of underground energy storage systems. The basic theory and 
much of the previous work in this area are found in the extensive 
review articles by Cheng [1] and Bejan [2]. Hany of these studies 
deal with steady or unsteady boundary layer flows along vertical 
and horizontal surfaces. 

The problem which deals with the flow of an incompressible, 
viscous fluid past an impulsively started infinite plate was first 
studied by Stokes [3] and the problem was later extended to bodies 
of different shapes as well as to viscous compressible fluids by a 
number of researchers. Closely related to the work reported in 
this chapter is the contribution of Soundalgekar [4], who has 
analysed the flow near an impulsively started infinite vertical 
plate, including the free convective effects caused by external 
heating or cooling of the plate. 

The main aim of the present analysis is to study the flow of 
a viscous fluid through a saturated porous medium adiacent to a 
moving vertical infinite plate. Considering the wall temperature 
to be different from that of the fluid, the resulting 



buoyancy-driven flow and the heat transfer in the fluid, are taken 
into account. The problem has an exact solution only for Pr ^ 1, 
Therefore, the solution has been obtained here by the approximate 
momentum- integral method for Pr < 1 and numerically by the finite 
difference method based on the Crank-Nicolson scheme. The results 
obtained by the two methods compare well. 

6.2 GOVERNING EQMATlOi^ 

Consider unsteady flow of an incompressible viscous fluid 
through a porous medium bounded by an infinite vertical moving 
wall. A co-ordinate system is introduced with x-axis taken along 
the wall in the upward direction and y-axis perpendicular to it. 
The flow through the porous medium is assumed to be governed by 
the generalized Darcy's law. It is also assumed that the 
Boussinesq approximation is applicable to the present problem. 
The governing equations for the problem are, 


5t 


sn 


CT - T^) t — 
&y 


(5.1) 

(5.2) 


where, u is the velocity in x-direction, t is the time, g is the 
acceleration due to gravity, /? is the coefficient of voluae 
expansion, K is the permeability of the porous medium, T is the 
temperature of the fluid, is the ambient temperature of the 
fluid, v is the kinematic viscosity and a is the thermal 
dif fusivity . 

The initial and boundary conditions are, 


t = 0 : u = 0, T = T^; y > 0, 
y = 0 : u = u^, T = T^; t > 0, 


(5.3) 



» fc.-t 


y-*<» : u-^0, T = T • t > 0, 

here, and are the tenperature and the velocity of the plate 
espect ively . 

The following non-dimensional variables are introduced. 


u 


Gr 



o 


vg/9 ) 

X — — (Grashof 


u 

o 



number) , 



(5.4) 


Pr = — (Prandtl number) and 
2 

11 

o 

Da = — j— K (Darcy number). 

V 

ibstituting (5.4) into (5.1) and (5.2), 
luations are obtained as. 


^U 

^t 

= + Gr & 

u 

" m 

SB 

1 


St 

~ U ' 



the non-dimensional 

(5.5) 

(5.6) 


le initial and boundary conditions reduce to , 
t = 0 : u = 0, d = 0, (y a 0), 

y = 0 : u = 1, e = 1, (t > 0), (5.7) 

y 00 : u 0 , d 0 , (t > 0). 

,3 METHC»:> OF SOLUTIOH 

Equation (5.6) has exact solution for all values of Pr, 
nereas the momentuiB equation (5.5) has exact solution for only Pr 

1. Using Laplace transfona technique, the solution of these 

:|uations is obtained [6] as follows: 


B = 




(5.8) 




u = -^ (1 - Gr Da) (exp [-y ] 


erfc (- 


24 t 




+ exp [y erfc + 4^)) + 


2^t 


-t/Da, 

r- ““ '-p^J 


« <«p t-Y j isj-f^irrTi cpfc - t -I > 

- .rfc c| - -iBr -ar-D ^ * Iv [)» ] 

* ■* tia (L-1) ■ ‘I " ■*Da 


+ Gr Da erfc (•^ ■J'^) . for Pr ^ 1 and 
u = -^ (1 - Gr Da) [exp [-y 4 erfc - ^ ^ ) 


(5-9) 


1 , a a y . . t 


m 


+ exp [y 4 u^3 erfc (-j^:^ + -I 


I5I 


)] 


+ Gr Da erfc ( 


TTt 


), for Pr = 1 


(5. JO) 


For Pr < 1, the equation (5.5) has no exact solution. 
Thus, the aain interest of the present study is to give the 
conplete solution of these equations. Consequently, the governing, 
equations are solved by finite difference aethod and the aoaentua 
integral aethod for Pr < 1 . 

Ca) Integral Method 

In this aethod, the velocity and teaperature profiles are 
chosen so as to vary saoothly froa the prescribed values at the 
surface and blend into the free streaa values at the edge of the 
boundary layer. A non-diaensional shape profile is assuaed and 
integrated over the boundary layer thickness, so that the local 
inaccuracies do not play a significant role. 

Let S, be the thickness of theraal boundary layer and S be 
the thickness of viscous boundary layer. 


Since the situation 



126 


under conalderatlon pertains to Pr < 1, it aay be concluded that 


<5 > S 

1 2 


Integrating equation (5.5) with respect to y, from y = 0 to 


, one gets , 


J dy = - c^)y=o J ^ J “ '^y- 


1 


.^u 


(5-11) 


0 0 0 
A third order profile is assumed for u in the form 

(5.12) 
in 

general, functions of t and are to be determined with the help of 


2 3 

u = b© * ^ 2 ^ * ^ 3 ^ ’ 

where, 17 = y /<5 The coefficients b , b. , b^ and b. are, 
Z 0 14 3 


boundary and initial conditions (5.7) and the smoothness condition 
= 0 at y = In addition to this, we use the following 

condition. 


^ . or « - 5 ^= 0 . 

d»y 


(5.13) 


which is obtained by evaluating the momentum equation at y = 0. 
The expression for u, thus takes the final form, 

A<5. 


3 ^13. 

„ = (1 _ 


- (17 - 217^ + T7^) , 


(5.14) 


where , 


A = 


1 - GrDa 
Da 


For a third order profile is assumed which can be expressed 


as : 


■■ 1 3 y 1 y- 

^ 2 <5, 2,3 


(5.15) 


where, <5^ = 2-^2 -/t/Pr. 

The above profile satisfies the temperature conditions of 

jlj lQ 

equation (5.7), the smoothness condition ^ = 0 at y « 6^ and 
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additional condition — * = 0, at y = 0 which reaulta fron 

ay 

evaluating equation (5.6) at y = 0. 

Subatituting (5.14) and (5.15) in (5.11) and integrating, one 

get 8: 


^ (18 - 3A ^ (72 + 12 A <5^^) + 36 Gr ^2t/Pr 

2 

^2 2 

~ Ua ^ ^2 ^ * (5.16) 

Thia equation can be reduced to the form, 

= {(144+24AZ) + 101.82337 Gr 4 (36-2A2 ) } / ( 18-3Az) , 

(5-17) 

where, z = ^ 


The values of z are obtained by the fourth-order Runge Kutta 
method. The expression for skin-friction, t, at the wall is given 


by: 


T 


-(— ) 
^ay^y=0 



(5.18) 


Cfc>3 Finite Difference Method 

To obtain the aolution by finite difference aethod, a region 
adjacent to the wall which ia aufficiently large to include both 
the viscoua and the thermal boundary layers is taken as the 
solution domain. Thia domain is divided into steps of equal size 
Ay in the y-direction, with grid points placed at the wall as well 
as the far stream boundary. In order to obtain the transient 
variation in the solution, time ia Incremented in equal steps of 
At. Since the dependent variables u and & are functions of two 
independent variables y and t, two subscripts are used to specify 
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the value of u and d at a given point at any tine, aay for 
instance, u(x , t.) = u. 

The difference equations corresponding to (5.5) and (5.6) 


are. 


At 


1 ''i-i,j ^''i,j “i+i,j 

_ I 

(Ay)^ 


(Ay) 


T 


+ Gr d. . - ■:rr=— [u. . + u. ...3, 

1 , J 2Da ‘ 1 , j 1, j + 1 


(5.19) 


d. - d. . , - 2d. . + e.^. . 

1,J + 1 l,J 1 f. 1-1, J i,J 1 + 1, J 

At " 2P? ^ 2 


^i-1, j+1 


(Ay) 
- 2d 


i,j+l * ®i+l,j+l 


]. (5.20) 


(Ay)‘ 


These equations lead to the following relat ionships between the 
values of u and B at six points (in space-tisie doiaain). 


At 

“X- u, . . . + (1 + 2X. + u, . “ X. u. . . . 

1 1-1, j+1 ^ I 2Da^ i,J+l 1 1+1, j+1 

= X, U. , . + (1-2X, - ^) u. . + X, u.^, . 

1 1-1, .) 1 2Da i,j 1 1 + 1 , J 


+ Gr d. .At, 
1 , J 


(5.21) 


-X., d. ^ . ^ + (1 + 2X.) d. - X-<9.^, ._^- 

2 1 - 1 , j +1 ^ 2 ' 1 , 1+1 2 1 + 1 , 1+1 


where , 


At ^ _ At 


(5.22) 


2(Ay) 


2 ’ 


^2 = 


2Pr (Ay)‘ 


The initial and boundary conditions 


1,0 

= 0 

and 

d. 

1,0 

= 0; 

for 

O, j 

= 1 

and 

d 

O, 1 

= 1; 

for 


take the following forn, 
all i, 

all j 0. (5.23) 
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, = 0 and B , * 0; for all j. 

) n, j 

where M la the grid point on the far atream boundary . 

At any time level, the above difference equationa will be 
aolved once for each point 1 i M-1, resulting in a system of 
M“1 simultaneoua equationa in M-1 unknowns, for each of the 
variables u and B. The equations exhibit a tridiagonal structure 
as follows: 

(1 + 2X, + u, . . - X, u, 

^ 1 2Da 1 , j+1 1 2 , j+1 


= "o.j + 1 ^ ■ 2^^ "l.j 

" “2.j ^ ^.j' 


(5.24) 


- S "i.j * N "i«.j 


+ Gr At d. 


for 2 i i lf-2 and 


(5.25) 


■^1 '"M-2,j + l + 2133^) “n-l.j + l 


1 M,j + 1 1 M-2,j 1 2Da n-l,j 

+ X- , + Gr At B . 

1 M,j M-l,j 


( 5 . 26 ) 


(1+2X2) 6 ^^ X2 ® 2 , j +1 


= X, d . + X, © . + (1-2X ) & + X 6 (5.27) 

2 o, J+1 2o,j 2 l,j 22, j 

. + fl + 2X_) e. ... - X^ 


-X B , . ^ (l + 2X_)d. 

2 i-l,j + l 2*' 1,1 + 1 2 i-i-l,j + l 


(5.28) 


'^2 ®M-2,j+i ( 1 + 2 X 2 ) - ^2 ^n,j+i 

^ ^2 V2,j ^ ^'-2^2^ Vl.J ^ ^2 ^H,j 


(5.29) 



The above system of equations have been solved by applying 
the tridiagonal matrix algorithm, once for the d-equations and 
then for the u-equations. 

5.4 RESULTS AND DISCUSSION 

For the problem under consideration, the solution domain 
extends from y = 0 to oo, in reality. However, calculations can be 
restricted to a smaller domain by observing that the analytical 
solutions of equations (5.5) and (5.6) which exist for Pr i 1 are 
close tozero around y .-x 4.0, and at all points after this. 
Similar results are also found in most of the flow problems over 
infinite plates. Reference in this regard may be made to the 
works of Soundalgekar [4] and Georgantopoulas et al . [5}. 

Therefore, although any value of y greater than 4.0 can be chosen 
to represent «>, in order to save computational time, it would be 
better to take an acceptably small value of y at which both u and 
& reduce to zero. In the present calculations, y = 5.0 is 

considered as the far-stream boundary. 

In the computations. Ay and At have been taken as 0.1 and 
0.005, respectively. Transient calculations have been continued 
till t = 1.0. The parameter values considered are: 

Pr = 0.71, 1-0, 7.0, Da = 0.25, 1.0, 4.0 and Gr = 1.0, 2.0. 

To check the correctness of taking y_^^= 5.0 as infinity, the 
program was run for y =6.0 and 7.0 and very little changes 

MBlX 

were observed in the temperature and velocity profiles, for the 
time values considered in the problem. To judge the accuracy of 
the results predicted, different values of Ay and At have been 
considered and it was verified that the step sizes employed in our 
calculations produce sufficiently accurate results. The 
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convergenc* of th« numerical scheme is thus tested. Finally, the 
skin ~ friction at the vail has been calculated by numerical 
differentiation, using Newton’s interpolation formula. 

Figure 1, shows the results obtained by the momentum integral 
method and the finite difference method using Crank-Nicolson 
scheme. The flow velocity along the plate has been plotted with 
respect to t, at differnt y-locations. The results obtained by 
the two methods compare well. It is also observed that the 

velocity tends to reach a steady value with time, although the 
time taken for such a stage increases with y. 

Rapt is and Singh [6] have obtained the velocity field using 

explicit differencing. Since it is well known that the explicit 

differencing in finite difference method produces only results 

which have an accuracy of 0(At), we have used the Crank-Nicolson 

scheme which is semi-implicit in nature. The Crank-Nicolson 

2 

scheme produces results which are 0(At ) accurate. There are, 
however, time step restrictions for both procedures to achieve 
numerically stable results. 

Figure 2, shows the velocity field, for different values of 
Da and Gr. Positive values of Gr have been considered which 
correspond to the case of heating of the plate. As expected, the 
velocity of the fluid increases with Gr, due to the greater role 
played by buoyancy forces. The velocity of the fluid increases 
with an increase in the value of permeability parameter also, 
since pore space increases with Da offering less resistance to 
flow. 

In figure 3, the variation in the velocity field for various 
values of Pr are shown. It is clear from the graph that the 



velocity at a point ia increased as the value of Pr is decreased. 
This may be attributed to the fact that at a lower value of Pr, 
the penetration of thermal effects is deeper due to larger thermal 
diffusivity. The increased thermal penetration results in an 
enhanced role for the buoyancy forces. As a consequence, the 
velocity values are larger. 

In figure 4, the transient velocity variation is shown for 
various Prandtl numbers. The value of velocity for every Pr starts 
from a minimum value and then increases with time as the fluid 
adjacent to the wall gets heated up. 

Table 1: Numerical values of skin-friction for Gr = 1.0, Pr = 0.71 


t 

Da 

T (Ref. [6]) 

T (Finite 
dif f erence) 

T (integral 
laethod) 

0.2 

1.0 

1.2447 

1.2360 

1.1654 

0.2 

4.0 

1.0596 

1.0828 

0.9574 

0.4 

1.0 

0.8664 

0.8631 

0.8117 

0.4 

4.0 

0.6030 

0.6331 

0.5171 


The values of skin-friction (dimensionless shear stress) obtained 
by various methods are compared in Table 1, for Gr = 1.0, Pr = 
0.71 and for different values of t and Da. The skin-friction 
obtained by the integral method differs from finite difference 
values by about 6% for Da = 1.0 while this difference becomes 
about 10% for Da = 4.0. 

It is observed that the value of t increases as Da decreases. 
This is due to the fact that the resistance to flow offered by the 
porous medium is higher at lower values of permeability, which in 


turn, results 

in 

steeper 

boundary 

layers 

leading to higher 

skin-friction. 

It 

is also 

observed 

that 

the skin-friction 
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decreases with tine, since the boundary layer thickness increases. 
G. S COKILUSION 

The free convection effects are studied on flow due to the 
inpulsive notion of an infinite vertical plate enbedded in a 
porous nediun. The problen has exact solution for Pr S: 1 , but for 
Pr < 1 , there is no exact solution. Here, results for Pr < 1 have 
been predicted by the approxinate nonentun- integral nethod and in 
general, by the finite difference nethod. The integral nethod 
solution is seen to be sufficiently accurate for low value of 
Darcy nunber. 
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Figure t Comparison of velocity profile 
by two methods. 




Figure 2 Velocity Profiles for Pr.0.71and for various values 
Of Da. 





Figure 3 Velocity Profiles for various 
values of Pr. 
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